Program VALUEFUNCTION

USE PARAMETERS
USE Functions_libs
USE My_ToolS
USE Interpolate_func
use omp_lib
Implicit NONE


! STATE VARIABLES

!     1. 18 for (educ) x (latent productivity type)* (latent health type) = "Type1" 
!     2. 8 for du (2) x dp (2) x s (2) = 8   = "Type2"
!     3. H (3) x R (2) = "Type 3"       
!     4. assets (12 points) 
!     5. human capital (experience) (5 grid points)
!     6. age
!     7. Marital Status (2)

! TYPE1 DEFINITION: educ x latent productivity x latent health - 16 possibilitis
!jtype1= 1 for e=1, zfix=1  , ht=1
!jtype1= 2 for e=1, zfix=1  , ht=2
!jtype1= 3 for e=1, zfix=2  , ht=1
!jtype1= 4 for e=1, zfix=2  , ht=2
!jtype1= 5 for e=1, zfix=3  , ht=1
!jtype1= 6 for e=1, zfix=3  , ht=2
!jtype1= 7 for e=2, zfix=1  , ht=1
!jtype1= 8 for e=2, zfix=1  , ht=2
!jtype1= 9 for e=2, zfix=2  , ht=1
!jtype1= 10 for e=2, zfix=2 , ht=2
!jtype1= 11 for e=2, zfix=3 , ht=1
!jtype1= 12 for e=2, zfix=3 , ht=2
!jtype1= 13 for e=3, zfix=1 , ht=1
!jtype1= 14 for e=3, zfix=1 , ht=2
!jtype1= 15 for e=3, zfix=2 , ht=1
!jtype1= 16 for e=3, zfix=2 , ht=2
!jtype1= 17 for e=3, zfix=3 , ht=1
!jtype1= 18 for e=3, zfix=3 , ht=2

! TYPE2 DEFINITION: Health shocks
!jtype2= 1 for du=1 , dp=1 , s=1 
!jtype2= 2 for du=1 , dp=2 , s=1 
!jtype2= 3 for du=1 , dp=1 , s=2 
!jtype2= 4 for du=1 , dp=2 , s=2  
!jtype2= 5  for du=2 , dp=1 , s=1 
!jtype2= 6  for du=2 , dp=2 , s=1  
!jtype2= 7  for du=2 , dp=1 , s=2 
!jtype2= 8  for du=2 , dp=2 , s=2 

! TYPE3 DEFINITION: combine H and R 
!jtype3=1 for jH=1, jR=1
!jtype3=2 for jH=1, jR=2
!jtype3=3 for jH=2, jR=1
!jtype3=4 for jH=2, jR=2
!jtype3=5 for jH=3, jR=1
!jtype3=6 for jH=3, jR=2

! TYPE4 DEFINITION: marital status and spousal work combined
!jtype4=1 for jmar=1  -- not married
!jtype4=2 for jmar=2, jos=1  -- married, spouse not working
!jtype4=3 for jmar=2, jos=2  -- married, spouse working

! Notes on notation
! jas - assets       ; jas_p - asssets next period
! je - education     ; jh - H
! jr - R             
! jx - human capital ; jfix - latent productivity type
! js,jdp,jdu  - all the health shocks      
! jht - latent health type 

Integer :: jo_lagged, joption_pay, jfix_par, jht_par, jpay, jhc_ret, jmc, jtrpay, jtype4, jtype4_p, jl_a, jm_a, ju_a,  jl_hc, jm_hc, ju_hc, jhc, jexp, jDI, jDI_p, jos, jos_p, jMEcat_p, jMEcat, jExperimentNumber, jspouse, counter, max_hk_age , jzhc_p,jzhc, jmar, jmar_p, jinc, jtype3,jtype3_p, mhk,jzt,jzt_p, kk, jas,jas_p,jsav,je,jh,jh_p, jr,jr_p,ji, jage,j, jht, jyears, jfix, jtype1,jtype2,jtype2_p,js,jdp,jdu,jx,jo,jo_p,  jins, jhrs , jx_age, jsav_upper_ret  , jsav_upper    , jinv_upper  , TR_indic , jhk


Real (Kind=dd), dimension (n_as,n_type1,n_type3,n_hc_grid,n_ret-1,n_o,n_mar,2)::  EXP_L_prime  !value functions; the last element is 1=not employed last period. 2=employed last period
Real (Kind=dd), dimension (n_as,n_type3,n_hc_grid,n_ret-1,n_o,n_mar,2)::  EXP_L_prime_simple   !value functions; the last element is 1=not employed last period. 2=employed last period 
Real (Kind=dd), dimension (n_as,n_type3,n_hc_grid,n_ret-1,n_zt,n_o,n_type4)::  EXP_L_prime_fix  !value functions  

Real (Kind=dd), dimension (n_o,n_type4,n_zt,n_as,n_type3,n_type1,n_hc_grid,n_ret-1)::  B !value function while working. 

Real (Kind=dd), dimension (n_o,n_MEcat,n_type4,n_zt,n_as,n_type3,n_type2,n_hc_grid,n_ret-1)::  B_1 !value function while working. (1=not have option to treat and not pay, 2=have all options)
Real (Kind=dd), dimension (n_o,n_MEcat,n_type4,n_zt,n_as,n_type3,n_type2,n_hc_grid,n_ret-1)::  B_2 !value function while working. (1=not have option to treat and not pay, 2=have all options)

Real (Kind=dd), dimension (n_trpay,n_o,n_MEcat,n_type4,n_zt,n_as,n_type3,n_type2,n_hc_grid,n_ret-1)::  B_fix !this is like B, but for one of the 18 fixed types. we want to work with smaller array so it's faster. at the end, combine this into B.

Real (Kind=dd), dimension (5,n_e,n_ht,n_type3,n_type2,n_MEcat,n_type2  ,n_type3)   :: prob_s_treat_new_64   
Real (Kind=dd), dimension (n_ret:n_age-1, n_e,n_ht,n_type3,n_type2, n_mar, n_MEcat, n_type2,  n_type3, n_mar )   :: prob_s_treat_new_ret

Real (Kind=dd)::  EXP_L_prime_INT, VF, V_death
Real (Kind=dd), dimension (2:3) :: VF_tr
Real (Kind=dd), dimension (n_sav,n_h,n_r,n_o,n_mar,n_zhc):: SAVE_EXP_L_prime_INT

Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_mar):: RET !value functions when retired
Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_ret : n_age,n_mar):: RET_PRIME !value functions when retired
Real (Kind=dd), dimension (n_MEcat,n_h,n_r,n_type2,n_mar):: RET_PRIME_INT

Real (Kind=dd), dimension (n_sav,n_MEcat,n_h,n_r,n_type2,n_mar):: SAVE_RET_PRIME_INT

Real  (Kind=dd) :: Resources, prob_s, prob_d, income_bt, income, base_1, base_2 , Total_tax, actual_earn, earnings_interpol, tax_SS, tax_med, earnings_capacity
Real (Kind=dd) :: last_value, maxB, valret   , check
Real (Kind=dd) :: UTI,x, EXPECTED_RET_PRIME , EXPECTED_B_PRIME_C_FLOOR_NT

Real (Kind=dd)::  slope_hc  , inter_hc    , slope_as, inter_as, HK_temp, CON, TR, ASP , MAX_SAV  , as_opt, TR_opt ,  INT_sav_upper_ret   , INT_sav_upper, EXPECTED_B_PRIME, EXPECTED_B_PRIME_C_FLOOR
Real (Kind=dd), dimension (n_o) ::  EXPECTED_B, jx_p_o
Real (Kind=dd), dimension (n_zhc,n_o) ::   HK_INTERPOL
Real (Kind=dd), dimension (n_e) :: scale_param

Real (Kind=dd), dimension (2) ::  income_pay, Total_tax_pay, base_2_pay, base_1_pay, CON_pay  ! these depend on whether you pay oop (1), or you don't pay oop (2)
Real (Kind=dd), dimension (2) :: MAX_SAV_pay, INT_sav_upper_pay  ! will calculate max savings for when you pay and don't pay OOP
Real (Kind=dd), dimension (2:3) :: prob_s_treat  ! this is used to calculate the probability of the state next period - 2 if you get treatement, 3 if you do not treat.
Real (Kind=dd), dimension (3) :: UTI_trpay, maxB_trpay,  EXPECTED_B_PRIME_trpay, last_value_trpay, V_cf_trpay, V_c_all_trpay
Real (Kind=dd), dimension (n_sav,3) :: SAVE_EXPECTED_B_PRIME_trpay
Integer, dimension (2) :: jsav_upper_pay

Real (Kind=dd), dimension (7) :: x_r
Real (Kind=dd), dimension (14) :: x_h   !     Contents of x: H_t0_Poor(1), H_t0_Fair(2), Dp(3), Du(4), TreatMC(1), TreatMC(2), TreatMC(3), TreatMC(4), Aget0(6), Aget0SQ(7), Aget0CB(8), EDU_COLG(9), EDU_SMCG(10)
REAL (Kind=dd) :: index_h
Real (Kind=dd), dimension (7) :: x_du
Real (Kind=dd), dimension (9) :: x_dp
Real (Kind=dd), dimension (7) :: x_s
Real (Kind=dd), dimension (n_e,n_hc_ret) :: x_ss

Integer, dimension (n_as,n_type1, n_ret, n_o, n_sav) :: AS_LB
Integer, dimension (n_hc_grid, n_e, n_type2, n_h, n_zhc, n_o, n_ret-1) ::  HK_LB  

integer m,n
real(dbl) yt(4), mx1, mx2
 
! Some useful variables and arrays for coding
! sav_grid: useful for interpolation to find the maximum savings
sav_grid(:)=0.0d0
 do jsav=1, n_sav
 sav_grid(jsav) = dble(jsav)
 end do

Call cpu_time(time1)
do jExperimentNumber = 0, 0  ! benchmark
    
 
!************************************************************************************************************************************************************ 
!       READING PARAMETERS, ESTIMATED FROM DATA, ASSUMED, AND INITIAL GUESSES
!************************************************************************************************************************************************************    

! Reading fixed assumed or estimated parameters: r, tax parameters, premiums
prem_ephi(:,:)=0.0d0 
prem_mcare(:)=0.0d0

open(202, file='Internal_Parameters/parameters_assumed.txt')
read (202,*) interest
read (202,*) taxc
read (202,*) tax_lambda(1) ! singles
read (202,*) tax_tau(1)
read (202,*) tax_lambda(2) ! married
read (202,*) tax_tau(2)
read (202,*) tau_ss
read (202,*) tau_med
read (202,*) yhat
read (202,*) penalty
read (202,*) prem_ephi(3,1)  ! singles premium
read (202,*) prem_ephi(3,3)  ! couples premium
read (202,*) prem_mcare(1) ! premium for a single person
close (202)


!*********  DEMOGRAPHICS AND SPOUSE *********!


! EQUIVALENCE SCALE BY AGE
equiv(:,:)=0.0d0
open(278, file='External_Parameters/parameters_fam_size.txt')
do je=1, n_e
do jage=1, n_ret-1 
read (278,*) equiv(je,jage)
equiv(je,jage) = 1.5d0 + (equiv(je,jage) -2.0d0 )* 0.3d0 ! Oxford scale. 
end do 
end do
  close (278)

 do jage= n_ret, n_age
   equiv(1,jage) = 1.5d0 
   equiv(2,jage) = 1.5d0 
   equiv(3,jage) = 1.5d0 
 end do
 
 
  
! EDUCATION
 p_e(:)=0.0d0
open(208, file='External_Parameters/Educ_dist.txt')
do je=1, n_e 
read (208,*) p_e(je)
end do 
close (208)


! SURVIVAL PROBABILITY
p_surv(:,:,:, :)=0.0 ! (n_h,n_e,n_mar,n_age)
open(195, file='External_Parameters/parameters_survival.txt')

 do jh=1, n_h ; do je=2,n_e; do jmar =1,n_mar ; do jage=26, n_age ! from age 50 (26 in the model)
   read (195,*) p_surv(jh,je,jmar,  jage)
 end do   ; end do ; end do   ; end do
 close (195)
 
 
 ! educations 1 and 2 are the same
   p_surv(:,1,:,  :)=p_surv(:,2,:,  :)

 
! fill in the rest of ages 1 to 25 ! assume linear trend starting from 0
   do jage=2, 25 ! prob is 0 at 25, so we have 24 positive numbers to add. so we divide the prob at jage=26 by 25 to get 25 equal intervals
   p_surv(:,:,:,  jage) = p_surv(:,:,:,  jage-1) + p_surv(:,:,:,  26)/ 25.0d0
   end do 
   ! change this to survival rather than mortality
   p_surv(:,:,:, :)=1.0d0-p_surv(:,:,:, :) 

   
! MARITAL STATUS, INITIAL PROBABILITIES AND TRANSITIONS   
! initial distribution of marital status at age 25
    In_mar(:, :)=0.0d0
open(193, file='External_Parameters/Marriage0.txt')
do je=1, n_e; do jh=1,n_h 
   read (193,*) In_mar(je, jh)
end do ; end do
close (193)

! MARITAL TRANSITIONS 
    ! M',M,je,jage,jh,jinc,jo
    Tp_mar(:,:,:,:,:,:,:)=0.0d0

! single to married
open(199, file='External_Parameters/Marriage1.txt')
do jage=1,n_ret-1;   do jo=1,n_o ;     do jinc=1, n_inc 
 read (199,*)  Tp_mar(2,1,1,jage,1,jinc,jo)
end do; end do; end do
close (199)
! e does not matter; h does not matter; 
Tp_mar(2,1,2,:,1,:,:) =Tp_mar(2,1,1,:,1,:,:)
Tp_mar(2,1,n_e,:,1,:,:) =Tp_mar(2,1,1,:,1,:,:)

Tp_mar(2,1,:,:,2,:,:) =Tp_mar(2,1,:,:,1,:,:)
Tp_mar(2,1,:,:,n_h,:,:) =Tp_mar(2,1,:,:,1,:,:)

! married to single
open(200, file='External_Parameters/Marriage2.txt')
 do jage=1,n_ret-1 
     do jh=1,n_h
         do jo=1,n_o
             do jinc=1, n_inc
                 do je=1,n_e 
 read (200,*)  Tp_mar(1,2,je,jage,jh,jinc,jo)
                 end do
             end do
         end do
     end do
 end do
close (200)
 
! retired individuals - no probability of marriage
Tp_mar(1,1,:,n_ret:n_age,:,:,:) =1.0d0
Tp_mar(2,1,:,n_ret:n_age,:,:,:) =0.0d0

! married to single - retired -- depends on je, jage only 
open(155, file='External_Parameters/Marriage3.txt')
do je=1,n_e 
do jage=n_ret, n_age-1
 read (155,*)  Tp_mar(1,2,je,jage,1,1,1)
 Tp_mar(1,2,je,jage,:,:,:)=Tp_mar(1,2,je,jage,1,1,1) 
end do
 end do
close (155)

Tp_mar(1,1,:,:,:,:,:)=1.0d0-Tp_mar(2,1,:,:,:,:,:)
Tp_mar(2,2,:,:,:,:,:)=1.0d0-Tp_mar(1,2,:,:,:,:,:) 




! PROBABILITY SPOUSE EMPLOYED
!  (n_ret-1,n_e,n_h,n_os) :: prob_os; ! spouse's employment: 1=not work, 2=work 
prob_os(:,:,:,:,:)=0.0d0

open(158, file='External_Parameters/Prob_spouse_work.txt')
 do jage=1,n_ret-1
     do je=1,n_e 
         do jh=1,n_h
             do jfix=1,n_fix
 read (158,*) prob_os(jage,je,jh,2,jfix) ! the last argument is productivity type.
            end do
         end do         
     end do
 end do
close (158)

prob_os(:,:,:,1,:) = 1.0d0 - prob_os(:,:,:,2,:) ! prob spouse does not work

! SPOUSE INCOME 

inc_spouse(:,:,:,:,:)=0.0d0 !(n_age,n_e,n_h,n_type4,n_fix) 
! spouse income when jtype4 is 1 is zero.
! at working ages, we allow for the spouse to have income even if she doesn't work. so if jtype=2, i.e., she does not work, she still has some low income. 

open(156, file='External_Parameters/Income_spouse1.txt')
 do jage=1,n_ret-1
     do je=1,n_e 
         do jh=1,n_h
             do jtype4=2,n_type4
                 do jfix=1,n_fix
 read (156,*)  inc_spouse(jage,je,jh,jtype4,jfix) ! last argument is jfix
               end do 
             end do
         end do         
     end do
 end do
close (156)



! for retired individuals, spouse's income only depends on education
! the spouse always has some income, doesn't depend on her work anymore, we just average over all spouses, so jtype4= 2 is the same as =3. 
! but for those not married, we keep at jtype4=1, that it's 0. we will use jmar in retirement, so when that's 1, income is 0. 


open(157, file='External_Parameters/Income_spouse_ret.txt')
     do je=1,n_e 
 read (157,*)  inc_spouse(n_ret,je,1,2,1)   ! jh=1, and jtype4=2   
 inc_spouse(n_ret+1:n_age,je,1,2,1)= inc_spouse(n_ret,je,1,2,1) ! same for all ages until n_age
 inc_spouse(n_ret:n_age,je,1,3,1)=inc_spouse(n_ret:n_age,je,1,2,1) ! same for type4=3 which is irrelevant in retirement since spouses don't work
 inc_spouse(n_ret:n_age,je,2,:,1)=inc_spouse(n_ret:n_age,je,1,:,1) ! doesn't depend on jh
 inc_spouse(n_ret:n_age,je,3,:,1)=inc_spouse(n_ret:n_age,je,1,:,1) ! doesn't depend on jh
     end do
close (157)

inc_spouse(n_ret:n_age,:,:,:,2) = inc_spouse(n_ret:n_age,:,:,:,1)
inc_spouse(n_ret:n_age,:,:,:,3) = inc_spouse(n_ret:n_age,:,:,:,1)

! SPOUSE MEDICAL OOP
!oop_spouse(jage,je,jtype4,jo)  
oop_spouse(:,:,:,:)=0.0d0  ! 0 will stay for those with jtype4=1 who are not married
open(159, file='External_Parameters/OOP_spouse.txt')
     do je=1,n_e 
         do jage=1, n_ret-1
read (159,*)  oop_spouse(jage,je,3,1) ! has insurance -- own or husband or both
read (159,*)  oop_spouse(jage,je,2,1) ! no insurance
oop_spouse(jage,je,2,3)=oop_spouse(jage,je,3,1) ! insured
oop_spouse(jage,je,2,5)=oop_spouse(jage,je,3,1) ! insured
oop_spouse(jage,je,3,2:5)=oop_spouse(jage,je,3,1) ! insured

oop_spouse(jage,je,2,2) = oop_spouse(jage,je,2,1) ! no insurance
oop_spouse(jage,je,2,4) = oop_spouse(jage,je,2,1) ! no insurance
         end do 
     end do
close (159)


open(169, file='External_Parameters/OOP_spouse_old.txt')
! the file is ages up to and including 83. in model years, that is 83-24=59
         do jage=n_ret, 59
read (169,*)  oop_spouse(jage,1,2,1) 
oop_spouse(jage,2:3,2,1)=oop_spouse(jage,1,2,1) ! All educ groups the same
         end do 
close (169)

! for ages 84+ , we just set oop equal to those at age 83. set all educ groups equal
oop_spouse(60:n_age,1,2,1)=oop_spouse(59,1,2,1)
oop_spouse(60:n_age,2,2,1)=oop_spouse(60:n_age,1,2,1)
oop_spouse(60:n_age,3,2,1)=oop_spouse(60:n_age,1,2,1)

! Income group thresholds
INC_TS(:)=0.0d0
open (232 , file='External_Parameters/parameters_inc_ts.txt')
read (232,*) INC_TS(1)
read (232,*) INC_TS(2)
read (232,*) INC_TS(3)
read (232,*) INC_TS(4)
close (232) 
 


! READ GUESSED VALUES
PT_penalty(:)=0.0d0
stdze(:)=0.0d0
varzfix(:)=0.0d0
varzt(:)=0.0d0
p_zhc(:,:,:)=0.0d0 
  param_calib(:)=0.0d0
FT_fc(:,:)=0.d0   ! disutility of work 
eta(:,:)=0.0d0
beta_hat(:)=0.0d0
theta_beq=0.0d0
k_beq=0.0d0
sig_beq=0.0d0
cons_floor(:,:,:,:)=0.0d0
c_floor_adjust_mar(:)=0.0d0
c_floor_age_adj_young(:)=0.0d0
c_floor_age_adj_old(:)=0.0d0
In_ht(:,:)=0.0d0
delta_o(:,:)=0.0d0
delta=0.0d0
param_o_trans(:,:)=0.0d0
leis_mar(:,:,:)=0.0d0
gamma=0.0d0
zfix(:,:)=0.0d0
   
open(203, file='Internal_parameters/parameters_internal_preferences.txt')
read (203,*) alpha
read (203,*) sigma
read (203,*) cost_death
do je=1, n_e
read (203,*) beta_hat(je)
end do
 do je=1, n_e; do jh= 1, n_h
 read (203,*)  FT_fc( je, jh)   ! disutility of work 
 end do ; end do
 read (203,*) eta(2,1) ! this is the stigma cost (disutility) of getting treated but not paying. Now this depends on insurance so we made a function of jo.
 read (203,*) eta(2,2)
 read (203,*) eta(2,3) ! has insurance
 read (203,*) eta(2,4)
 read (203,*) eta(2,5) ! has insurance
do jo=1,5,2 ; do je=1,3      
  read (203,*)    leis_mar(je,jo,1) ! adjustment to time endowment if married and wife does not work
end do ; end do 
do jo=1,5,2 ; do je=1,3      
  read (203,*)    leis_mar(je,jo,2) ! adjustment to time endowment if married and wife works
end do ; end do 
read (203,*) theta_beq
read (203,*) k_beq
read (203,*) sig_beq
close (203)

open(204, file='Internal_parameters/parameters_internal.txt')
do je=1, n_e
read (204,*)	  PT_penalty(je)
read (204,*)      stdze(je)
read (204,*)      varzfix(je)
read (204,*)      varzt(je)
read (204,*)	  b0(je)	
read (204,*)	  b1(je)	
read (204,*)	  b2(je)
read (204,*)	  b3(je)	
read (204,*)	  b4(je)
read (204,*)	  b5(je)
end do
do je=1, n_e
    read (204,*)	  b6(je)
end do

do je=1, n_e
read (204,*) cons_floor(1,1,je, 1) ! Poor H; first element is age=1 
read (204,*) cons_floor(1,1,je, 2) ! Not Poor H
end do
do je=1, n_e
read (204,*) c_floor_adjust_mar(je) ! we need a higher consumption floor for the married, even after we adjust for equivalence scale. this parameter is a scaling parameter.
end do
do je=1, n_e
read (204,*) c_floor_age_adj_young(je)
read (204,*) c_floor_age_adj_old(je)
end do
do je=1, n_e ; do jfix=1, n_fix
read (204,*) In_ht(je,jfix)  
end do; end do
do je=1, n_e 
read (204,*) param_o_trans(je,1) 
read (204,*) param_o_trans(je,2) 
read (204,*) param_o_trans(je,3) 
read (204,*) param_o_trans(je,4) 
read (204,*) param_o_trans(je,5) 
read (204,*) param_o_trans(je,6) 
end do

 read (204,*) delta_o(1,1) 
  read (204,*) delta_o(2,1)
   read (204,*) delta_o(3,1)
    read (204,*) delta_o(1,2) 
     read (204,*) delta_o(2,2)
      read (204,*) delta_o(3,2)
read (204,*) zfix(1,1)
read (204,*) zfix(1,2) 
read (204,*) zfix(1,n_fix) 
read (204,*) zfix(2,1)
read (204,*) zfix(2,2)
read (204,*) zfix(2,n_fix) 
read (204,*) zfix(3,1)
read (204,*) zfix(3,2)
read (204,*) zfix(3,n_fix) 
close (204)


! same leisure cost, ESHI does not matter.
do je=1,3 ; do jos=1,2
leis_mar(je,2,jos)=leis_mar(je,3,jos)
leis_mar(je,4,jos)=leis_mar(je,5,jos)
end do; end do


!Fill in HK shocks
p_zhc(:,3,:) = p_zhc(:,2,:) ! equal to prob for being employed for all o's where employed.
p_zhc(:,4,:) = p_zhc(:,2,:)
p_zhc(:,5,:) = p_zhc(:,2,:)

! Consumption floor
!singles - second element =1 for singles
do je=1, n_e
cons_floor(1:n_age,1,je, 1) =cons_floor(1,1,je, 1) !DI
cons_floor(1:n_age,1,je, 2) =cons_floor(1,1,je, 2) !fair H
cons_floor(1:n_age,1,je, 3) =cons_floor(1,1,je, 2) ! cons floor for those in good H. not getting DI
end do

!married - depends on age because of equivalence scale (second element =2)
do je=1, n_e
do jage=1, n_ret-1
cons_floor(jage,2,je, 1) = cons_floor(jage,1,je, 1)* equiv(je,jage)  ! POOR H - DO NOT SCALE UP ANY FURTHER FOR MARRIED
cons_floor(jage,2,je, 2) = cons_floor(jage,1,je, 2)* equiv(je,jage)*(c_floor_adjust_mar(je)) ! FAIR H 
cons_floor(jage,2,je, 3) = cons_floor(jage,1,je, 3)* equiv(je,jage)*(c_floor_adjust_mar(je)) ! GOOD H
    end do 
end do 

do je=1, n_e
do jage=n_ret, n_age
    do jh=1, n_h
cons_floor(jage,2,je, jh) = cons_floor(jage,1,je, jh)* equiv(je,jage)
end do
    end do 
end do 

! adjusting the consumption floor by age
do jage=21, n_ret-1; do je=1, n_e
   cons_floor(jage,2,je, 2:3) = cons_floor(jage,2,je, 2:3)+   cons_floor(jage,2,je, 2:3)* (dble(jage)-20.0d0)*c_floor_age_adj_old(je)
end do ; end do 


do jage=19, 1, -1; do je=1, n_e
   cons_floor(jage,2,je, 2:3) = cons_floor(jage,2,je, 2:3)-   cons_floor(jage,2,je, 2:3)* (20.0d0 - dble(jage))*c_floor_age_adj_young(je)
end do ; end do


do je=1, n_e; do jh=1, n_h; do jmar=1,n_mar; do jage=1, n_age
cons_floor_at(jage,jmar,je,jh)    = cons_floor(jage,jmar,je,jh)*(1.0d0+taxc) ! consumption floor including tax
end do; end do; end do; end do


! Hours of work offered
hrs_offer(:)=0.0d0
hrs_offer(1)=0.0d0
hrs_offer(2)=hrs_pt  
hrs_offer(3)=hrs_pt 
hrs_offer(4)=hrs_ft 
hrs_offer(5)=hrs_ft 

! Wage penalty for working PT
!convenient to have this as function of o
wage_penalty(:,:)=0.0d0
wage_penalty(1,:) = 0.0d0
wage_penalty(2,1) = PT_penalty(1)
wage_penalty(2,2) = PT_penalty(2)
wage_penalty(2,3) = PT_penalty(3)
wage_penalty(3,1) = PT_penalty(1)
wage_penalty(3,2) = PT_penalty(2)
wage_penalty(3,3) = PT_penalty(3)
wage_penalty(4,:) = 1.0d0
wage_penalty(5,:) = 1.0d0

do je=1, n_e
b1(je)=b1(je)/2.d0
b2(je)=b2(je)/4.d0
b3(je)=b3(je)/8.d0
end do

! Shocks to human capital
zhc(:)=0.0d0
zhc(1)=1.0d0



! HUMAN CAPITAL GRID
HK(:)=0.0d0
do jx=0, n_hc ! from 0 to 39
HK(jx) = dble(jx)
end do

! determine max possible human capital at every age
max_age_HC(:)=0.0d0
max_age_HC(1) = 0.0d0
max_age_HC(2) = zhc(n_zhc)  

do jage=3, n_ret-1
max_age_HC(jage) = (max_age_HC(jage-1) +1.0d0) *zhc(n_zhc)  ! maximum yrs of HC at every age
end do 

!now determine the grid
!xfn - this function gives the actual number of years of experience given any grid point from 1 to 5 of human capital for every age
xfn(:, :)=0.0d0


do jage=1, n_ret-1
xfn(1, jage)=0.d0 ! first grid point always corresponds to 0 human capital.
xfn(2, jage)=  max_age_HC(jage) * 0.22d0  !  
xfn(3, jage)=  max_age_HC(jage) * 0.45d0    ! 
xfn(4, jage)=  max_age_HC(jage) * 0.65d0    ! 
xfn(5, jage)= max_age_HC(jage)
end do



! Social Security Income 
pen(:,:,:)=0.0d0 


open(187, file='Internal_parameters/SS_inc.txt')
do je =1, n_e; do jfix=1, n_fix; do jhc_ret =1, n_hc_ret
read (187,*)	  pen(jhc_ret, je, jfix)
end do ; end do; end do 
close(187)

!Category thresholds for human capital at age 65
hc_TS(:,:)=0.0d0
open(68, file='Internal_Parameters/parameters_hc_ts.txt')
do je =1, n_e
 read (68,*) hc_TS(1,je)
  read (68,*) hc_TS(2,je)
   !read (68,*) hc_TS(3)
  end do
close (68)

! TYPE1 distribution (fixed types) AND LATENT SKILL
p_fix(:,:)=1.0d0/3.0d0 ! probability of a latent skill type within an education gruop is 1/3. We have 3 skill groups in each.

In_type1(:)=0.0d0
In_type1(1)=p_e(1)*(1.0d0/3.0d0)*In_ht(1,1)
In_type1(2)=p_e(1)*(1.0d0/3.0d0)*(1.0d0-In_ht(1,1))
In_type1(3)=p_e(1)*(1.0d0/3.0d0)*In_ht(1,2)
In_type1(4)=p_e(1)*(1.0d0/3.0d0)*(1.0d0-In_ht(1,2))
In_type1(5)=p_e(1)*(1.0d0/3.0d0)*In_ht(1,3)
In_type1(6)=p_e(1)*(1.0d0/3.0d0)*(1.0d0-In_ht(1,3))

In_type1(7)=p_e(2)*(1.0d0/3.0d0)*In_ht(2,1)
In_type1(8)=p_e(2)*(1.0d0/3.0d0)*(1.0d0-In_ht(2,1))
In_type1(9)=p_e(2)*(1.0d0/3.0d0)*In_ht(2,2)
In_type1(10)=p_e(2)*(1.0d0/3.0d0)*(1.0d0-In_ht(2,2))
In_type1(11)=p_e(2)*(1.0d0/3.0d0)*In_ht(2,3)
In_type1(12)=p_e(2)*(1.0d0/3.0d0)*(1.0d0-In_ht(2,3))

In_type1(13)=p_e(3)*(1.0d0/3.0d0)*In_ht(3,1)
In_type1(14)=p_e(3)*(1.0d0/3.0d0)*(1.0d0-In_ht(3,1))
In_type1(15)=p_e(3)*(1.0d0/3.0d0)*In_ht(3,2)
In_type1(16)=p_e(3)*(1.0d0/3.0d0)*(1.0d0-In_ht(3,2))
In_type1(17)=p_e(3)*(1.0d0/3.0d0)*In_ht(3,3)
In_type1(18)=p_e(3)*(1.0d0/3.0d0)*(1.0d0-In_ht(3,3))

!EMPLOYMENT OFFERS PARAMETERS ********

! Probability of the employment offer at age 25 by educ.
Pr_o_25(:,:)=0.0d0
open(105, file='Internal_Parameters/parameters_offers_25.txt')
do je=1, n_e; do jo=1, n_o
   read (105,*) Pr_o_25(jo, je)  
end do; end do 
close (105)

Pr_o_i(:,:,:,:)=0.0d0
do je=1, n_e !  jo_p, jo, je, jage
Pr_o_i(2,1,je,1) = (1.0d0 - param_o_trans(je,3)  ) * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(2,2,je,1) = (1.0d0 - param_o_trans(je,2)  ) * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(2,3,je,1) = (1.0d0 - param_o_trans(je,2)  ) * (1.0d0 - param_o_trans(je,4)  )
Pr_o_i(2,4,je,1) = (1.0d0 - param_o_trans(je,1)  ) * (1.0d0 - param_o_trans(je,6)  )
Pr_o_i(2,5,je,1) = (1.0d0 - param_o_trans(je,1)  ) * (1.0d0 - param_o_trans(je,4)  )

Pr_o_i(3,1,je,1) = (1.0d0 - param_o_trans(je,3)  ) *  param_o_trans(je,5)  
Pr_o_i(3,2,je,1) = (1.0d0 - param_o_trans(je,2)  ) *  param_o_trans(je,5)  
Pr_o_i(3,3,je,1) = (1.0d0 - param_o_trans(je,2)  ) *  param_o_trans(je,4)  
Pr_o_i(3,4,je,1) = (1.0d0 - param_o_trans(je,1)  ) *  param_o_trans(je,6)  
Pr_o_i(3,5,je,1) = (1.0d0 - param_o_trans(je,1)  ) *  param_o_trans(je,4)  

Pr_o_i(4,1,je,1) =  param_o_trans(je,3)   * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(4,2,je,1) =  param_o_trans(je,2)   * (1.0d0 - param_o_trans(je,5)  )
Pr_o_i(4,3,je,1) =  param_o_trans(je,2)   * (1.0d0 - param_o_trans(je,4)  )
Pr_o_i(4,4,je,1) =  param_o_trans(je,1)   * (1.0d0 - param_o_trans(je,6)  )
Pr_o_i(4,5,je,1) =  param_o_trans(je,1)   * (1.0d0 - param_o_trans(je,4)  )

Pr_o_i(5,1,je,1) =  param_o_trans(je,3)   *  param_o_trans(je,5)  
Pr_o_i(5,2,je,1) =  param_o_trans(je,2)   *  param_o_trans(je,5)  
Pr_o_i(5,3,je,1) =  param_o_trans(je,2)   *  param_o_trans(je,4)  
Pr_o_i(5,4,je,1) =  param_o_trans(je,1)   *  param_o_trans(je,6)  
Pr_o_i(5,5,je,1) =  param_o_trans(je,1)   *  param_o_trans(je,4)  

end do 

! ages 25-53 are the same
do jage=2,n_ret-1
Pr_o_i(:,:,:,jage) =   Pr_o_i(:,:,:,1)
end do


! HIGH SCHOOL
! now allow this to depend on age. people 54 and over have a positive probability of having no offer, and this increases linearly in age 
scale_param(:)=0.0d0
je=1
do jage=30, 35
    do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-29))*delta_o(je,1) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(1,jo,je,jage)=min(1.0d0,dble((jage-29))*delta_o(je,1) + Pr_o_i(1,jo,je,jage) ) ! about 4% per year
! reduce all other probabilities proportionally
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage)* scale_param(je))
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* scale_param(je))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
end do
end do 

! from age 60 onwards, we have a steeper increase in the probability of no offer
do jage=36, n_ret-1
     do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-29))*delta_o(je,2) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(1,jo,je,jage)=min(1.0d0,dble((jage-29))*delta_o(je,2) + Pr_o_i(1,jo,je,jage) ) ! about 4% per year
! reduce all other probabilities proportionally
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage)* scale_param(je))
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* scale_param(je))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
end do
end do



! SOME COLLEGE
! allow this to depend on age. people 54 and over have a positive probability of having no offer, and this increases linearly in age 
scale_param(:)=0.0d0
je=2
do jage=23, 32 ! 
    do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-22))*delta_o(je,1) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .85d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .85d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
    end do
end do 

! from age 58 onwards, we have a steeper increase in the probability of no offer
do jage=33, n_ret-1
     do jo=1,5
scale_param(je)=  (1.0d0 - dble((jage-23))*delta_o(je,2) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .9d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .9d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
end do
end do

! COLLEGE
! allow this to depend on age. people 54 and over have a positive probability of having no offer, and this increases linearly in age 
scale_param(:)=0.0d0
je=3
do jage=23, 32 
    do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-22))*delta_o(je,1) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))

Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .85d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .85d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
end do 
end do

! from age 60 onwards, we have a steeper increase in the probability of no offer
do jage=33, n_ret-1 ! 57
     do jo=1,5
    scale_param(je)=  (1.0d0 - dble((jage-23))*delta_o(je,2) - Pr_o_i(1,jo,je,jage)) /(1.0d0 - Pr_o_i(1,jo,je,jage))
Pr_o_i(4,jo, je,jage)= max(0.0d0, Pr_o_i(4,jo, je,jage)* scale_param(je))
Pr_o_i(5,jo, je,jage)= max(0.0d0, Pr_o_i(5,jo, je,jage)* scale_param(je))
Pr_o_i(2,jo, je,jage)= max(0.0d0, Pr_o_i(2,jo, je,jage) * .9d0**(dble((jage-22)))  )
Pr_o_i(3,jo, je,jage)= max(0.0d0, Pr_o_i(3,jo, je,jage)* .9d0**(dble((jage-22)))  )
Pr_o_i(1,jo,je,jage) = 1.0d0 -Pr_o_i(2,jo, je,jage)-Pr_o_i(3,jo, je,jage)-Pr_o_i(4,jo, je,jage)-Pr_o_i(5,jo, je,jage)
end do
end do


!Premiums
prem_mcare(2) =prem_mcare(1) * 2.0d0 ! jmar =1 if single and 2 if married
prem_ephi(5,3)=prem_ephi(3,3) ! couples premium 

! single or spouse not insured
prem_ephi(3,2)=prem_ephi(3,1) ! singles
prem_ephi(5,1)=prem_ephi(3,1) ! singles
prem_ephi(5,2)=prem_ephi(3,1) ! singles

! only spouse has insurance but not husband
prem_ephi(1,3)=prem_ephi(3,1)
prem_ephi(2,3)=prem_ephi(3,1)
prem_ephi(4,3)=prem_ephi(3,1)

! ******************************** HEALTH PROCESS ******************************
! DISTRIBUTION OF H, R AT AGE 25
In_h(:, :,:)=0.0d0
In_r(:, :)=0.0d0


open(201, file='Internal_Parameters/parameters_H_25.txt')
do je=1, n_e; do jht=2,1,-1; do jh=1,n_h 
   read (201,*) In_h(je, jh,jht)
end do ; end do; end do
close (201)



open(207, file='External_Parameters/parameters_R25.txt')
   read (207,*) In_r(1, 2) ! reading probabilities of high risk
   In_r(1, 1)=1.0d0 - In_r(1, 2) 
   In_r(2, :)=In_r(1, :) 
   In_r(3, :)=In_r(1, :)
close (207)


! H Transitions and Probability shocks and R not observed
 param_h_cutoff_A(:)=0.0d0
 param_h_cutoff_B(:)=0.0d0
 param_h(:)=0.0d0
 neu_shocks(:,:,:)=1.0d0 ! if you are treated, then prob =1. we only change this later if not treated 
 neu_R=1.0d0
 Shock_type(:,:,:,:)=0.0d0
 treat_age_adj=0.0d0
 open(22, file='Internal_parameters/parameters_internal_health.txt')
  

! Parameters for H
   read (22,*) param_h_cutoff_A(1) ! Cutoff parameters for Type A -- good type
   read (22,*) param_h_cutoff_A(2)
   read (22,*) param_h_cutoff_B(1) ! Cutoff parameters for Type B
   read (22,*) param_h_cutoff_B(2)
read (22,*) param_h(1)   ! H_t0_Poor(1), H_t0_Fair(2) 
read (22,*) param_h(2)
read (22,*) param_h(3) ! Dp(3), Du(4), TreatMC(5), Aget0(6), Aget0SQ(7), Aget0CB(8), EDU_COLG(9), EDU_SMCG(10)
read (22,*) param_h(4)
read (22,*) param_h(5) ! for MC
read (22,*) param_h(6) ! for MC
read (22,*) param_h(7) ! for MC
read (22,*) param_h(8) ! for MC
read (22,*) param_h(9)
read (22,*) param_h(10)
read (22,*) param_h(11)
read (22,*) param_h(12)
read (22,*) param_h(13)
read (22,*) param_h(14) ! indicator for age 65+
read (22,*) neu_shocks(1,1,1) !
read (22,*) neu_shocks(2,1,1) 
read (22,*) neu_shocks(3,1,1) 
read (22,*) neu_shocks(1,2,1) 
read (22,*) neu_shocks(2,2,1) 
read (22,*) neu_shocks(3,2,1) 
read (22,*) neu_shocks(1,3,1) 
read (22,*) neu_shocks(2,3,1) 
read (22,*) neu_shocks(3,3,1) 
read (22,*) neu_R(1) 
read (22,*) Shock_type(1,3,1,1) !insured so jo=3 or 5 ! cannot treat, insured
read (22,*) Shock_type(1,3,2,1) !insured so jo=3 or 5  ! cannot treat, insured
read (22,*) Shock_type(1,3,3,1) !insured so jo=3 or 5  ! cannot treat, insured
read (22,*) Shock_type(1,1,1,1) ! uninsured so jo=1,2,4 ! cannot treat, uninsured
read (22,*) Shock_type(1,1,2,1) ! uninsured so jo=1,2,4  ! cannot treat, uninsured
read (22,*) Shock_type(1,1,3,1) ! uninsured so jo=1,2,4  ! cannot treat, uninsured
read (22,*) Shock_type(2,3,1,1) ! can treat, but no option to default
read (22,*) Shock_type(2,3,2,1) ! can treat, but no option to default
read (22,*) Shock_type(2,3,3,1) ! can treat, but no option to default
read (22,*) Shock_type(2,1,1,1) ! can treat, but no option to default
read (22,*) Shock_type(2,1,2,1) ! can treat, but no option to default
read (22,*) Shock_type(2,1,3,1) ! can treat, but no option to default
read (22,*) treat_age_adj ! parameter to lower probability of "can't treat" at older ages
close (22)

! fill in for all jo in Shock_type
! uninsured
Shock_type(:,2,:,1) = Shock_type(:,1,:,1)
Shock_type(:,4,:,1) = Shock_type(:,1,:,1)
! insured
Shock_type(:,5,:,1) = Shock_type(:,3,:,1)
! Fix shock_type for the last category where you can treat and not pay. 
Shock_type(3,:,:,1) = 1.0d0 - Shock_type(1,:,:,1) -Shock_type(2,:,:,1)
! fix the ages
! from 25-44, they stay as they are, but from age 45 (21 in model) the probability that you cannot treat decreases linearly. the other two go up proportionally
do jage=1,20
Shock_type(:,:,:,jage) = Shock_type(:,:,:,1) 
end do 

do jage=21, n_ret-1
    Shock_type(:,:,:,jage) = Shock_type(:,:,:,1)
end do 

!we only adjust if uninsured and not in poor H - lower probabilities of "cannot treat" at old ages. 
do jage=21, n_ret-1 ; do jh=2,3
Shock_type(1,1,jh,jage) = Shock_type(1,1,jh,1)- treat_age_adj*(dble(jage)-20.0d0) ! cannto treat
Shock_type(1,2,jh,jage) = Shock_type(1,2,jh,1)- treat_age_adj*(dble(jage)-20.0d0) ! cannto treat
Shock_type(1,4,jh,jage) = Shock_type(1,4,jh,1)- treat_age_adj*(dble(jage)-20.0d0) ! cannto treat

Shock_type(2,1,jh,jage) =( Shock_type(2,1,jh,jage) / (Shock_type(2,1,jh,jage) + Shock_type(3,1,jh,1))) * (1.0d0 - Shock_type(1,1,jh,jage) ) ! now make these proportional 
Shock_type(2,2,jh,jage) =( Shock_type(2,2,jh,jage) / (Shock_type(2,2,jh,jage) + Shock_type(3,2,jh,1))) * (1.0d0 - Shock_type(1,2,jh,jage) )
Shock_type(2,4,jh,jage) =( Shock_type(2,4,jh,jage) / (Shock_type(2,4,jh,jage) + Shock_type(3,4,jh,1))) * (1.0d0 - Shock_type(1,4,jh,jage) )
end do ; end do

Shock_type(3,:,:,:) = 1.0d0 - Shock_type(1,:,:,:) -Shock_type(2,:,:,:) !  making sure all add up to 1 - this is for having all 3 options including treat and not pay

!Neu_R for uninsured
neu_R(2)=neu_R(1)
neu_R(4)=neu_R(1)

! H transitions
 x_h(:)=0.0d0  
 index_h=0.0d0
 
 Tp_h(:,:, :, :, :,:,:,:)=0.0d0
 trprob_h(:,:,:,:,:,:,:,:)=0.0d0
 
 ! IF TREATED
 
do jh = 1, 3 ; do je = 1, 3 ; do jage = 1,  n_ret-1 ; do jdp = 0, 1 ; do jdu = 0, 1 
!x_h   !     Contents of x: H_t0_Poor(1), H_t0_Fair(2), Dp(3), Du(4), MC=1(5), MC=2(6), MC=3(7), MC=4(8), Aget0(9), Aget0SQ(10), Aget0CB(11), EDU_COLG(12), EDU_SMCG(13)

      x_h = (/ 0.0d0, 0.0d0, DBLE(jdp), DBLE(jdu), 0.0d0, 0.0d0, 0.0d0, 0.0d0, DBLE(jage), DBLE((jage)**2), DBLE((jage)**3), 0.0d0, 0.0d0, 0.0d0 /)

    IF (jh==1) THEN
        x_h(1)=1.0d0
    ELSEIF (jh==2) THEN
        x_h(2)=1.0d0
    ENDIF
    IF (je==3) THEN
        x_h(12)=1.0d0
    ELSEIF (je==2) THEN
        x_h(13)=1.0d0
    ENDIF
    If (jage.ge.41) then
        x_h(14)=1.0d0
    end if 
    
    
        index_h=DOT_PRODUCT(param_h(:),x_h(:))

        trprob_h(jh,1,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1))) ! last argument is =1 (i.e., getting treated, or no shocks)
        trprob_h(jh,2,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        
        trprob_h(jh,1,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))
        trprob_h(jh,2,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))

end do; end do; end do; end do; end do


x_h(:) =0.0d0
index_h=0.0d0

! if not treated ! i'm only doing this for ages up to and including 64. the rest get "treated"
do jh = 1, 3 ; do je = 1, 3 ; do jage = 1, n_ret-1 ; do jdp = 0, 1 ; do jdu = 0, 1 ; do jmc = 2, 5
!x_h   !     Contents of x: H_t0_Poor(1), H_t0_Fair(2), Dp(3), Du(4), MC=1(5), MC=2(6), MC=3(7), MC=4(8), Aget0(9), Aget0SQ(10), Aget0CB(11), EDU_COLG(12), EDU_SMCG(13)
    x_h = (/ 0.0d0, 0.0d0, DBLE(jdp), DBLE(jdu), 0.0d0, 0.0d0, 0.0d0, 0.0d0, DBLE(jage), DBLE((jage)**2), DBLE((jage)**3), 0.0d0, 0.0d0 ,0.0d0 /) 

    IF (jh==1) THEN
        x_h(1)=1.0d0
    ELSEIF (jh==2) THEN
        x_h(2)=1.0d0
    ENDIF
    IF (je==3) THEN
        x_h(12)=1.0d0
    ELSEIF (je==2) THEN
        x_h(13)=1.0d0
    ENDIF
    if (jmc==2) then
       x_h(5)=1.0d0 
       else if (jmc==3) then
       x_h(6)=1.0d0 
        else if (jmc==4) then
       x_h(7)=1.0d0 
        else if (jmc==5) then
       x_h(8)=1.0d0 
    end if 
    
        index_h=DOT_PRODUCT(param_h(1:13),x_h(1:13))

        trprob_h(jh,1,je,1,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        trprob_h(jh,2,je,1,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        
        trprob_h(jh,1,je,2,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))
        trprob_h(jh,2,je,2,jage,jdp+1,jdu+1,jmc) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))

end do; end do; end do; end do; end do; end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! now shift the probabilities at 65+
param_h_cutoff_A(1)=param_h_cutoff_A(1)+ 0.6d0
param_h_cutoff_B(1)=param_h_cutoff_B(1)+ 0.6d0


x_h(:) =0.0d0
index_h=0.0d0
 ! IF TREATED (all treat after 65)
 
do jh = 1, 3 ; do je = 1, 3 ; do jage = n_ret,  n_age ; do jdp = 0, 1 ; do jdu = 0, 1 
!x_h   !     Contents of x: H_t0_Poor(1), H_t0_Fair(2), Dp(3), Du(4), MC=1(5), MC=2(6), MC=3(7), MC=4(8), Aget0(9), Aget0SQ(10), Aget0CB(11), EDU_COLG(12), EDU_SMCG(13)
      x_h = (/ 0.0d0, 0.0d0, DBLE(jdp), DBLE(jdu), 0.0d0, 0.0d0, 0.0d0, 0.0d0, DBLE(jage), DBLE((jage)**2), DBLE((jage)**3), 0.0d0, 0.0d0, 0.0d0 /)

    IF (jh==1) THEN
        x_h(1)=1.0d0
    ELSEIF (jh==2) THEN
        x_h(2)=1.0d0
    ENDIF
    IF (je==3) THEN
        x_h(12)=1.0d0
    ELSEIF (je==2) THEN
        x_h(13)=1.0d0
    ENDIF
    If (jage.ge.41) then
        x_h(14)=1.0d0
    end if 
    
    
        index_h=DOT_PRODUCT(param_h(:),x_h(:))

        trprob_h(jh,1,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1))) ! last argument is =1 (i.e., getting treated, or no shocks)
        trprob_h(jh,2,je,1,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_A(1)))
        
        trprob_h(jh,1,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))
        trprob_h(jh,2,je,2,jage,jdp+1,jdu+1,1) = 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(2))) - 1.0d0 / (1.0d0 + DEXP(index_h-param_h_cutoff_B(1)))

end do; end do; end do; end do; end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!






trprob_h(:,3,:,:,:,:,:,:) = 1.0d0 - trprob_h(:,1,:,:,:,:,:,:) - trprob_h(:,2,:,:,:,:,:,:)

do jh_p = 1,n_h ; do jh = 1,n_h ; do je=1,n_e ; do jht=1,2 ;  do jage = 1, n_age; do jmc=1,4 ! 2=treat, we will have 4 grid points for charges.they are all the same if get treated.
Tp_h(jh_p,jh,je,jht,jage,1,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,1)   ! no shock ! in trprob_h, jmc=1 if treat so this is why here, the last argument is always 1 because we're considering the case for treated
Tp_h(jh_p,jh,je,jht,jage,3,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,1)   ! only s

Tp_h(jh_p,jh,je,jht,jage,6,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,1) ! both shocks du and dp
Tp_h(jh_p,jh,je,jht,jage,8,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,1) ! both shocks du and dp

Tp_h(jh_p,jh,je,jht,jage,2,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,1) ! dp only
Tp_h(jh_p,jh,je,jht,jage,4,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,1) ! dp only

Tp_h(jh_p,jh,je,jht,jage,5,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,1) ! du only
Tp_h(jh_p,jh,je,jht,jage,7,2,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,1) ! du only

end do; end do; end do; end do; end do; end do

! those over 65+ are "treated, so fill in the array for them
Tp_h(:,:,:,:,n_ret:n_age,:,1,:) = Tp_h(:,:,:,:,n_ret:n_age,:,2,:)


! if not treated - only possible up to age 64, inclusive.
do jh_p = 1,n_h ; do jh = 1,n_h ;do je=1,n_e ; do jht=1,2 ;  do jage = 1, n_ret-1 ; do jmc=1,4 ! 2=treat
Tp_h(jh_p,jh,je,jht,jage,1,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,1)   ! no shock ! here, the person doesn't get penalized for not getting treated.
Tp_h(jh_p,jh,je,jht,jage,3,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,1,jmc+1)   ! only s

Tp_h(jh_p,jh,je,jht,jage,6,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,jmc+1) ! both shocks du and dp
Tp_h(jh_p,jh,je,jht,jage,8,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,2,jmc+1) ! both shocks du and dp

Tp_h(jh_p,jh,je,jht,jage,2,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,jmc+1) ! dp only
Tp_h(jh_p,jh,je,jht,jage,4,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,2,1,jmc+1) ! dp only

Tp_h(jh_p,jh,je,jht,jage,5,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,jmc+1) ! du only
Tp_h(jh_p,jh,je,jht,jage,7,1,jmc) = trprob_h(jh,jh_p,je,jht,jage,1,2,jmc+1) ! du only

end do; end do; end do; end do; end do; end do


! write to file to plot
1009 Format(I1,1x,I1,1x,I1,1x,I1,1x,I2,1x,I1,1x,I1,1x,I1,1x, F9.5)
open(234, file='Simulated_Data/H_transitions.txt')
do jh_p = 1,n_h ; do jh = 1,n_h ; do je=1,n_e ; do jht=1,2 ;  do jage = 1, n_age ; do jmc=1,4 ; do jtype2=1,n_type2; do jpay=1,2
write(234,1009) jh_p, jh, je, jht, jage, jmc, jtype2 , jpay, Tp_h(jh_p,jh,je,jht,jage,jtype2,jpay,jmc)
end do; end do; end do; end do; end do; end do; end do; end do
close(234)

! du, dp, s
shocks_risk(:,:, :, :, :)  =0.d0 !(n_e,n_type2, n_h, n_r, n_age)

open(222, file='External_Parameters/parameters_health_shocks.txt')
do je=1,n_e; do jtype2=1, n_type2 ; do jh=1,n_h ; do jr=1,n_r;  do jage=1, n_age
   read (222,*) shocks_risk(je,jtype2,jh,jr,jage)
end do ; end do ; end do  ; end do ; end do
close (222) 



! R transitions
p_risk(:,:,:,:)=0.d0 !(n_r,n_r,n_age, n_h)
open(223, file='External_Parameters/parameters_health_R.txt')  
do jh=1,n_h ; do jr=1,n_r;  do jage=1, n_age
   read (223,*) p_risk(jr,2, jage,jh)
end do ; end do ; end do  
close (223) 

p_risk(:,1,:,:) = 1.0d0 - p_risk(:,2,:,:) 


! ****************** MEDICAL TREATMENT COSTS ******************
!OOP
oop(:, :, :,:,:)=0.0d0

! ages<65
! for those with ESHI insurance. 
! non-catastrophic
open(198, file='External_Parameters/ME_ins_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (198,*) oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+1,1,5)=oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+2,1,5)=oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+3,1,5)=oop(jtype2, jh, jage,1,5)
   oop(jtype2, jh, jage+4,1,5)=oop(jtype2, jh, jage,1,5)
end do ; end do ; end do  
close (198)   

oop(:, :, :,1,3) =oop(:, :,:,1,5) !  insured

! Insurance and catastrophic
open(168, file='External_Parameters/ME_ins_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (168,*) oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+1,2,5)=oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+2,2,5)=oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+3,2,5)=oop(jtype2, jh, jage,2,5)
   oop(jtype2, jh, jage+4,2,5)=oop(jtype2, jh, jage,2,5)
end do ; end do ; end do  
close (168)   

oop(:, :, :,2,3) =oop(:, :,:,2,5) !  insured

! age<65, no ESHI 
! non-catastrophic
open(45, file='External_Parameters/ME_no_ins_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (45,*) oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+1,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+2,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+3,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+4,1,1)=oop(jtype2, jh, jage,1,1)
end do ; end do ; end do  
close (45)   

oop(:, :, :,1,2) =oop(:, :,:,1,1) ! not insured
oop(:, :, :,1,4) =oop(:, :,:,1,1) ! not insured

! Insurance and catastrophic
open(46, file='External_Parameters/ME_no_ins_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (46,*) oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+1,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+2,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+3,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+4,2,1)=oop(jtype2, jh, jage,2,1)
end do ; end do ; end do  
close (46)   

oop(:, :, :,2,2) =oop(:, :,:,2,1) ! not insured
oop(:, :, :,2,4) =oop(:, :,:,2,1) ! not insured

! ages >=65
open(47, file='External_Parameters/ME_ret_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=41, 51,5
   read (47,*) oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+1,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+2,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+3,1,1)=oop(jtype2, jh, jage,1,1)
   oop(jtype2, jh, jage+4,1,1)=oop(jtype2, jh, jage,1,1)
end do ; end do ; end do  
close (47)  


open(48, file='External_Parameters/ME_ret_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=41, 51,5
   read (48,*) oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+1,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+2,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+3,2,1)=oop(jtype2, jh, jage,2,1)
   oop(jtype2, jh, jage+4,2,1)=oop(jtype2, jh, jage,2,1)
end do ; end do ; end do  
close (48) 

! fill in ages 80 and over
do jage=56,n_age,1
oop(:, :, jage,:,1)=oop(:, :, 55,:,1)
end do


! for retired - jo doesn't matter 
oop(:, :, n_ret:n_age,:,2) = oop(:, :, n_ret:n_age,:,1)
oop(:, :, n_ret:n_age,:,3) = oop(:, :, n_ret:n_age,:,1)
oop(:, :, n_ret:n_age,:,4) = oop(:, :, n_ret:n_age,:,1)
oop(:, :, n_ret:n_age,:,5) = oop(:, :, n_ret:n_age,:,1)

! Medical Charges 

! these are never used in the code. we are just writing them to file in the end. 
! the purpose is simply in terms of matching moments based on Charges rather than OOP. 
! these moments are needed only for those <65

! dimension(n_type2, n_h, n_ret-1, n_MEcat) :: Charges

! non-catastrophic
open(49, file='External_Parameters/Charges_1.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (49,*) Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+1,1)=Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+2,1)=Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+3,1)=Charges(jtype2, jh, jage,1)
   Charges(jtype2, jh, jage+4,1)=Charges(jtype2, jh, jage,1)
end do ; end do ; end do  
close (49) 

! catastrophic
open(50, file='External_Parameters/Charges_2.txt')
do jtype2=1, n_type2 ; do jh=1,n_h ;  do jage=1, 36,5
   read (50,*) Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+1,2)=Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+2,2)=Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+3,2)=Charges(jtype2, jh, jage,2)
   Charges(jtype2, jh, jage+4,2)=Charges(jtype2, jh, jage,2)
end do ; end do ; end do  
close (50) 

!Categories for Medical Charges - read thresholds
open(67, file='External_Parameters/parameters_mc_ts.txt')
 read (67,*) MC_TS(1)
  read (67,*) MC_TS(2)
   read (67,*) MC_TS(3)
close (67)

! PROBABILITY OF CATASTROPHIC ME 
Prob_MEcat(:)=0.0d0
Prob_MEcat(1) = 0.95d0   
Prob_MEcat(2) = 1.0d0- Prob_MEcat(1) 


!*************** SICK DAYS    ****************

phi(:,:,:)=0.0d0  ! (educ, H, type2)  
 ! those without any shocks have 0 sick days
! sick days are the same for Some College and College. They are also the same for Fair and Poor H. 
open(230, file='External_Parameters/sick_days.txt')
do jtype2=1, n_type2
read (230,*) phi(1,1,jtype2) ! HS, Fair/Poor H
phi(1,2,jtype2)=phi(1,1,jtype2)
end do 

do jtype2=1, n_type2
read (230,*) phi(1,3,jtype2) ! HS, Good H
end do 

do jtype2=1, n_type2
read (230,*) phi(2,1,jtype2) ! S Col, Fair/Poor H
phi(2,2,jtype2)=phi(2,1,jtype2) ! Fair H the same
phi(3,1,jtype2)=phi(2,1,jtype2) ! College is the same as S College
phi(3,2,jtype2)=phi(2,2,jtype2) ! College is the same as S College
end do 

do jtype2=1, n_type2
read (230,*) phi(2,3,jtype2) ! S Col, Good H
phi(3,3,jtype2)=phi(2,3,jtype2) ! College the same
end do 

close (230)

! Sick days, or leisure cost of illness for non-workers
phi_nw(:,:)=0.0d0
open(233, file='External_Parameters/sick_days_nw.txt')
do jh=2,n_h
do jtype2=1, n_type2
read (233,*) phi_nw(jh,jtype2) 
end do 
end do 
close (233)
phi_nw(1,:) =phi_nw(2,:) ! we combined poor and fair health

! Define leisure conditional on o, educ and health 
! account for the additional disutility of work.; this takes into account sick days (lost leisure) for non-workers. Note that for workers, the sick days are included in work hours - they just don't get paid for them. 
leis(:,:,:,:)=0.0d0
do je=1, n_e ; do jh=1, n_h; do jtype2=1, n_type2
leis(je,1,jh,jtype2) = 1.0d0 - phi_nw(jh,jtype2) ! not working get full leisure
leis(je,2:3,jh,jtype2)=  1.0d0 - hrs_pt - .75*FT_fc( je, jh)  ! we are subtracting the same fixed cost for both PT and FT. this is because we need to discourage PT from working. 
leis(je,4:5,jh,jtype2)=  1.0d0 - hrs_ft - FT_fc( je, jh) !FT
end do ; end do ; end do


!**************** WAGES AND BASE EARNINGS ****************
stdzt(:)=0.0d0
stdzt(1)=varzt(1)**0.5d0
stdzt(2)=varzt(2)**0.5d0
stdzt(3)=varzt(3)**0.5d0


Call Deterministic_z ! This subroutine gives zd

zt(:,:)=0.0d0
do je=1, n_e
zt(je,1)=-1.d0*stdzt(je) 
zt(je,n_zt)=1.d0*stdzt(je) 
Call discretize_normal(zt(je,1),zt(je,n_zt),n_zt,0.d0,stdzt(je),zt(je,:),p_zt(je,:,2)) 
Call discretize_normal(zt(je,1),zt(je,n_zt),n_zt,b6(je),stdzt(je),zt(je,:),p_zt(je,:,1)) 

end do

 z(:,:,:,:,:)=0.0d0
 base_earn(:,:,:,:,:,:)=0.0d0 
do je=1, n_e; do jfix=1, n_fix;  do jzt=1, n_zt 

do jx=0,n_hc;  do jh=1,n_h 
  z(jx,jh,je,jfix,jzt) = exp(zd(je,jx,jh)+zfix(je,jfix)+zt(je,jzt))  ! this gives wage offers
  
  base_earn(jx,jh,je,jfix,jzt,1)=0.0d0 ! not working
  base_earn(jx,jh,je,jfix,jzt,2:3)= z(jx,jh,je,jfix,jzt)*wage_penalty(2,je)     ! part time
  base_earn(jx,jh,je,jfix,jzt,4:5)= z(jx,jh,je,jfix,jzt)  ! FT
end do
end do; end do ; end do  ; end do



! Read initial assets
assets_initial(:)=0.0d0
as(:,:, :)=0.0d0 


as(1,:,:)=0.d0
as(2,:, :)=10000.d0/5200.d0
as(3,:, :)=20000.d0/5200.d0
as(4,:, :)=50000.d0/5200.d0
as(5,:, :)=75000.d0/5200.d0
as(6,:, :)=100000.d0/5200.d0 
as(7,:, :)=135000.d0/5200.d0  
as(8,:,:)=200000.d0/5200.d0 
as(9,:, :)=400000.d0/5200.d0 
as(10,:, :)=750000.d0/5200.d0
as(11,:, :)=1100000.d0/5200.d0
as(12,:, :)=1700000.d0/5200.d0

! the discrete savings choice
sav_ret(:)=0.0d0

open(231, file='External_Parameters/savings_grid_ret.txt')
do jsav=1, n_sav
   read (231,*) sav_ret(jsav)
end do 
close (231)
sav_ret(:)=sav_ret(:)/5200.0d0


sav(:,:)=0.0d0
open(252, file='External_Parameters/savings_grid.txt')
do jsav=1, n_sav
   read (252,*) sav(5,jsav)
end do 
close (252)
sav(5,:)=sav(5,:)/5200.0d0


! those not working and those without insurance have to be able to dissave a lot more bc of possibility of catastrophic ME
open(253, file='External_Parameters/savings_grid_PT.txt')
do jsav=1, n_sav
   read (253,*) sav(1,jsav)
end do 
close (253)
sav(1,:)=sav(1,:)/5200.0d0
sav(2,:)=sav(1,:)
sav(3,:)=sav(1,:)
sav(4,:)=sav(1,:)

! FIND THE LOWER BOUND GRID POINT FOR ASSET INTERPOLATION

 ASP=0.0d0   
 AS_LB(:,:,:,:,:)=0.0d0
    do jas=1, n_as
        do jtype1=1, n_type1 
            do jage=1, n_ret-1
                do jo=1, n_o
                    do jsav=1,n_sav
       ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)  
        
        if (ASP.ge.as(n_as,jtype1,jage+1)) then 
            AS_LB(jas,jtype1,jage+1,jo,jsav)= n_as
        else
            jl_a=1
            ju_a=n_as
        do
            if (ju_a-jl_a <= 1) exit
            jm_a = int(( ju_a + jl_a ) / 2)
            if (ASP >= as(jm_a,jtype1,jage+1) ) then
	            jl_a = jm_a
            else
	            ju_a = jm_a
            end if
        end do
        end if 
        
       AS_LB(jas,jtype1,jage+1,jo,jsav)=jl_a  ! note I'm writing jage+1 here; this is for next age grid for when we interpolate for future period
                    end do
                end do
            end do 
        end do
    end do
    
! FIND THE LOWER BOUND GRID POINT FOR HUMAN CAPITAL INTERPOLATION

    !careful because this is the point on the grid tomorrow at t+1 , but we have written jage not jage+1. 
    
 HK_temp=0.0d0
 HK_LB(:,:,:,:,:,:,:)=0.0d0
 
  
    do jage=1, n_ret-2 
    do jhc=1, n_hc_grid ! if this will be constant and not depend on age, i.e., 5 points for all ages, replace this with that later.
     do je=1, n_e
         do jtype2=1, n_type2
            do jh=1, n_h
                 jx_p_o(:)=0.0d0
                 jx_p_o(1) = xfn(jhc,jage)
                 jx_p_o(2) =min( (xfn(jhc,jage)+( max((hrs_pt - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft)))  ,  dble(n_hc) ) 
                 jx_p_o(3) = jx_p_o(2)
                 jx_p_o(4) = min( (xfn(jhc,jage)+( max((hrs_ft - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft)))  ,  dble(n_hc) )
                 jx_p_o(5) = jx_p_o(4)
        
        do jzhc=1,n_zhc ! these are shocks that happen just before next period, so it's going to be jzhc_p
            do jo=1, n_o !current jo
          HK_temp  = max(0.0d0, min( jx_p_o(jo) *zhc(jzhc), dble(n_hc) )) ! n_hc is 0-100
          !!!!!!!!!!
         if (HK_temp.ge. xfn(n_hc_grid, jage+1) ) then 
            HK_LB(jhc,je,jtype2,jh,jzhc,jo,jage+1)= n_hc_grid  ! note again this is jage+1
        else
            jl_hc=1
            ju_hc=n_hc_grid
        do
            if (ju_hc-jl_hc <= 1) exit
            jm_hc = int(( ju_hc + jl_hc ) / 2)
            if (HK_temp >= xfn(jm_hc   , jage+1) ) then
	            jl_hc = jm_hc
            else
	            ju_hc = jm_hc
            end if
            HK_LB(jhc,je,jtype2,jh,jzhc,jo,jage+1)=jl_hc
        end do
        end if  

            end do     
        end do  

            end do 
         end do
     end do
    
    end do 
    end do
    
! Multiplication of probabilities  -  "prob_s_treat_new"
  prob_s_treat_new(:,:,:,:,:,:,:,:,:,:)=0.0d0  
  do jage=1,n_ret-1; do jo=1,n_o; do je=1,n_e; do jht=1, n_ht; do jh=1,n_h; do jr=1, n_r; do jtype2=1, n_type2 
       do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)
             call TYPE3_DET(jh,jr,jtype3)
                          
        
  !YOU GET TREATED 
        prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) = Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh) 
  ! Not treated      
        prob_s_treat_new(2,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,1)*p_risk(jr, jr_p,jage,jh) 
        prob_s_treat_new(3,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,2)*p_risk(jr, jr_p,jage,jh) 
        prob_s_treat_new(4,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,3)*p_risk(jr, jr_p,jage,jh) 
        prob_s_treat_new(5,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) =  Pr_o_i(jo_p,jo, je,jage)*Tp_h(jh_p, jh,je,jht,jage, jtype2,1,4)*p_risk(jr, jr_p,jage,jh) 
 
            
       end do; end do ; end do ; end do ; end do  ; end do ; end do
  end do; end do ; end do ; end do 
 
  !Probability at 64 (n_ret-1)
 prob_s_treat_new_64(:,:,:,:,:,:,:,:)=0.0d0    ! (treat and MC  , je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p)
    do je=1,n_e; do jht=1, n_ht; do jh=1,n_h; do jr=1, n_r; do jtype2=1, n_type2 
        do jh_p=1, n_h ; do jr_p=1, n_r ; do jtype2_p=1, n_type2 ; do jMEcat_p=1, n_MEcat 
                        call TYPE3_DET(jh_p,jr_p,jtype3_p)
             call TYPE3_DET(jh,jr,jtype3)
            
        
!treat
            prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,2,1)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
!not treat . 2= lowest MC category, and 5 highest.
            prob_s_treat_new_64(2,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,1)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
            prob_s_treat_new_64(3,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,2)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
            prob_s_treat_new_64(4,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,3)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)
            prob_s_treat_new_64(5,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) =   Tp_h(jh_p, jh,je,jht,n_ret-1, jtype2,1,4)*p_risk(jr, jr_p,n_ret-1,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, n_ret)*Prob_MEcat(jMEcat_p)

        
        end do; end do ; end do ; end do ; end do  ; end do ; end do; end do ; end do
        
        
        
 ! Probabilities at 65+
      prob_s_treat_new_ret(:,:,:,:,:,:,:,:,:,:)=0.0d0    ! (jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p)
 do jage=n_ret, n_age-1;   do je=1,n_e; do jht=1, n_ht; do jh=1,n_h; do jr=1, n_r; do jtype2=1, n_type2 ; do jmar=1,n_mar
        do jh_p=1, n_h ; do jr_p=1, n_r ; do jtype2_p=1, n_type2 ; do jMEcat_p=1, n_MEcat ; do jmar_p=1,n_mar
                        call TYPE3_DET(jh_p,jr_p,jtype3_p)
             call TYPE3_DET(jh,jr,jtype3)   
prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) = Tp_h(jh_p, jh,je,jht,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) * Tp_mar(jmar_p,jmar,je,jage,1,1,1) 
        end do; end do ; end do ; end do ; end do  ; end do ; end do
        end do ; end do   ; end do ; end do ; end do
        
!***************************************************************************
! Initialize value function - AGE 100
!***************************************************************************
    ! in the last period, you decide optimal bequests ! iterate to find optimal
jage=n_age  !dead for sure at n_age+1
RET_PRIME(:,:,:, :,:,:,:,:,:)   =0.d0

do jhc_ret = 1, n_hc_ret ! social security depends on earning capacity at age 64
do je = 1, n_e
do jfix=1,min(n_fix, I_state_fix)  ! we don't need this, but filling up array

do jht=1,min(n_ht, I_state_ht)  ! we don't need this, but filling up array
call TYPE1_DET(je,jht,jfix,jtype1)
do jas=1,n_as
do jMEcat=1, n_MEcat  
do jmar=1,n_mar
do jdu=1,n_du ; do jdp=1,n_dp ; do js=1,n_s 
call TYPE2_DET(jdu,jdp,js,jtype2)         
do jh=1, n_h ; do jr=1, n_r
    
! we write  inc_spouse(jage,je,1,jmar) because here jtype4 is 1 if not married and 2 if married ; jh does not matter so we just write 1.
Resources = interest*as(jas,jtype1,n_age) -prem_mcare(jmar) - oop(jtype2,jh,n_age,jMEcat,1) - oop_spouse(n_age,je,jmar,1) +pen(jhc_ret,  je,jfix) +  inc_spouse(n_age,je,1,jmar,jfix)  - max(0.0d0,  (( ((max(0.0d0, (interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret, je,jfix) +  inc_spouse(jage,je,1,jmar,jfix)-prem_mcare(jmar) - max(0.d0, oop(jtype2,jh,n_age,jMEcat,1) + oop_spouse(n_age,je,jmar,1)  - 0.075d0*((interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret,je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )))*5200.0d0) - tax_lambda(jmar) *( ((max(0.0d0, (interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret, je,jfix) +  inc_spouse(jage,je,1,jmar,jfix)-prem_mcare(jmar) - max(0.d0, oop(jtype2,jh,n_age,jMEcat,1) + oop_spouse(n_age,je,jmar,1)  - 0.075d0*((interest-1.0d0)*as(jas,jtype1,n_age) + pen(jhc_ret,je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )))*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 ))
CON=0.0d0
ASP=0.0d0
TR=0.0d0
V_cf_trpay(:)=0.0d0
V_c_all_trpay(:) =0.0d0

if (Resources.le.cons_floor_at(jage, jmar,je,jh)) then
   CON= cons_floor(jage, jmar,je,jh)  !this is the cons floor
   call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
  RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=UTI + beta_hat(je)* ( cost_death + Bequest_val(0.0d0) )
else 
    
! Calculate value of consuming the consumption floor and leaving everything else as bequest.
    call func_uti_ret(je,jh,jmar,jtype2,cons_floor(jage, jmar,je,jh),UTI)
V_cf_trpay(1) = UTI + beta_hat(je)* ( cost_death + Bequest_val(Resources-cons_floor_at(jage, jmar,je,jh)) )
! calculate value of consuming all and saving nothing
    call func_uti_ret(je,jh,jmar,jtype2,(Resources/(1.0d0+taxc)),UTI)
V_c_all_trpay(1) = UTI+ beta_hat(je)* ( cost_death + Bequest_val(0.0d0) )

    !find max savings
    MAX_SAV=0.0d0
    INT_sav_upper_ret=0.0d0
    jsav_upper_ret=0.0d0
    MAX_SAV=Resources-cons_floor_at(jage, jmar,je,jh) -as(jas,jtype1,n_age)
    
        
                            jsav_upper_ret=0
          do while ( MAX_SAV .ge.  sav_ret(jsav_upper_ret+1))
           jsav_upper_ret=jsav_upper_ret +1
           if (jsav_upper_ret.EQ.n_sav) exit
          end do 
          

               
        
! LOOP FOR SAVINGS ***************************
        maxB=V_cf_trpay(1)
    do jsav= jsav_upper_ret , 1, -1
    CON= ( Resources -  as(jas,jtype1,n_age) - sav_ret(jsav) )/(1.0d0+taxc)
    ASP=  as(jas,jtype1,n_age)+  sav_ret(jsav)
!Get values
    if (ASP.lt.0.d0)  exit
    if (CON.lt.cons_floor(jage, jmar,je,jh)) then 
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
    end if  
        call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
	    RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=    UTI+ beta_hat(je)* ( cost_death + Bequest_val(ASP) )
  
        
        if (RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar).lt.maxB) exit
			  	    maxB=RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)                  
   
    end do !for jsav for retired ppl
    
     RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=maxB 
     
  ! now find the maximum considering the extremes
    if (RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar).eq.0.0d0) then
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
    end if 
    
  
        RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar)=max(RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,n_age,jmar),  V_c_all_trpay(1)) 
   
         
end if ! for being on the cons floor

                      
end do; end do !For jr, jh
end do ! jmar
end do; end do; end do !js, jdp, jdu
end do ! jMEcat
end do !For jas
end do ! jht
end do ! jfix
end do !for  je
end do !jhc_ret

!***************************************************************************************************************************************
!  Value function - Retirees
!***************************************************************************************************************************************
   ! (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_mar):: RET


do jage=n_age-1,n_ret,-1  
   PRINT *,'Current Age:',jage
   RET(:,:,:,:,:,:,:,:)=0.d0   
jfix=1 
do je = 1, n_e
do jht=1, min(n_ht, I_state_ht)    ! if we want only one state, set I_state=1, and we do only jht=1
    
call TYPE1_DET(je,jht,jfix,jtype1)

do jdu=1,n_du ; do jdp=1,n_dp 
    js=1 ! this used to be "do js=1,n_s" but we don't need to loop here since s does not matter for next period variables. so leave the s loop for after. Do calcualtion only for s=1 and then fill in D_prime because s=2 will be the same as for s=1
    
call TYPE2_DET(jdu,jdp,js,jtype2) 
do jh=1, n_h ; do jr=1, n_r  
    
    call TYPE3_DET(jh,jr,jtype3)
    
do js=1,n_s 
    call TYPE2_DET(jdu,jdp,js,jtype2) ! now call this to get the rigth type2 given s=1 and s=2
do jas=1,n_as
do jMEcat=1, n_MEcat  
do jmar=1,n_mar 

do jhc_ret = 1, n_hc_ret ! social security depends on earning capacity at age 64

! Initiate    
   TR=0.d0
   CON=0.d0
   ASP=0.d0
   UTI=0.d0
   TR_indic=0
   MAX_SAV=0.0d0
   V_cf_trpay(:)=0.0d0
   V_c_all_trpay(:) =0.0d0
    
! HH income:
income = max(0.0d0, (interest-1.0d0)*as(jas,jtype1,jage) + pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -  prem_mcare(jmar)  - max(0.d0,  oop(jtype2, jh, jage,jMEcat,1) +  oop_spouse(jage,je,jmar,1) - 0.075d0*((interest-1.0d0)*as(jas,jtype1,jage) + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )) !taxable income. the last term is OOP in excess of a threshold of your income, which is deductable.
base_2 = prem_mcare(jmar) + oop(jtype2, jh, jage,jMEcat,1) + oop_spouse(jage,je,jmar,1)  + max(0.0d0,  (( (income*5200.0d0) - tax_lambda(jmar) *( (income*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) !! the last thing is the tax fn that was =Tax_inc(jmar,income) ! total resources subtracted: premiums, OOP, and tax, at the HH level
! we write  inc_spouse(jage,je,1,jmar) because here jtype4 is 1 if not married and 2 if married, so this works for retirees

!determine if HH gets transfer
if (((interest*as(jas,jtype1,jage)+ pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) .lt. (cons_floor_at(jage, jmar,je,jh)) ) then                                                                                                                                      
    TR= cons_floor_at(jage, jmar,je,jh) - (interest*as(jas,jtype1,jage)+ pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2 ) 
    CON=cons_floor(jage, jmar,je,jh)
    TR_indic=1 
    MAX_SAV=0.0d0 
else
   TR=0.d0
   CON=0.d0
   TR_indic=0
   MAX_SAV=  interest*as(jas,jtype1,jage)  + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) - base_2 -cons_floor_at(jage, jmar,je,jh) -as(jas,jtype1,jage)
end if 

!*****  Big IF statement for TR
if (TR_indic.eq.1) then
!solve for RET   
    RET_PRIME_INT(:,:,:,:,:)=0.0d0
    EXPECTED_RET_PRIME=0.d0 
    V_death=0.0d0
    V_death=cost_death +Bequest_val(0.0d0) 
    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
        call TYPE3_DET(jh_p,jr_p,jtype3_p)
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) *( p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  + (1.0d0- p_surv(jh_p,je,jmar_p, jage+1))  * V_death )
    end do ; end do  ; end do; end do; end do
 
        call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
        RET(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,jmar)=UTI+ beta_hat(je)*EXPECTED_RET_PRIME 
    
else if (TR_indic.eq.0) then
    
! Calculate value of consuming the consumption floor and saving all else: ASP= MAX_SAV+as(jas,jtype1,jage) 
    RET_PRIME_INT(:,:,:,:,:)=0.0d0
    EXPECTED_RET_PRIME=0.d0 
    ASP=0.0d0
    ASP=MAX_SAV+as(jas,jtype1,jage)
    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
        call TYPE3_DET(jh_p,jr_p,jtype3_p)
                call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) *( p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)  + (1.0d0- p_surv(jh_p,je,jmar_p, jage+1))  * (cost_death +Bequest_val(ASP) ) )
    end do ; end do  ; end do; end do; end do
    
    call func_uti_ret(je,jh,jmar,jtype2,cons_floor(jage, jmar,je,jh),UTI)
  V_cf_trpay(1) = UTI + beta_hat(je)* EXPECTED_RET_PRIME

! Calculate value of consuming all and saving nothing
  EXPECTED_RET_PRIME=0.0d0
      do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
          call TYPE3_DET(jh_p,jr_p,jtype3_p)
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) *( p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  + (1.0d0- p_surv(jh_p,je,jmar_p, jage+1))  * (cost_death +Bequest_val(0.0d0) ) )
    end do ; end do  ; end do; end do; end do
    call func_uti_ret(je,jh,jmar,jtype2,((((interest*as(jas,jtype1,jage)+ pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) )/(1.0d0+taxc)),UTI) ! cash on hand divided by cons tax
V_c_all_trpay(1) = UTI+ beta_hat(je)* EXPECTED_RET_PRIME


        
        
                     jsav_upper_ret=0
          do while ( MAX_SAV .ge.  sav_ret(jsav_upper_ret+1))
           jsav_upper_ret=jsav_upper_ret +1
           if (jsav_upper_ret.EQ.n_sav) exit
          end do 
          
  
               
               
               
           
! LOOP FOR SAVINGS 
        maxB=V_cf_trpay(1) ! initialize to max savings off the grid
    do jsav= jsav_upper_ret , 1, -1
!Solve value 
    TR=0.d0
    CON= ( ((interest-1.0d0))*as(jas,jtype1,jage) + inc_spouse(jage,je,1,jmar,jfix) + pen(jhc_ret,  je,jfix) -base_2  - sav_ret(jsav) )/(1.0d0+taxc)
    ASP=  as(jas,jtype1,jage)+  sav_ret(jsav)
    
    if (ASP.lt.0.d0)  exit
    if (CON.lt.cons_floor(jage, jmar,je,jh)) then 
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
    end if  
           
            EXPECTED_RET_PRIME=0.d0
            RET_PRIME_INT(:,:,:,:,:)=0.0d0
             do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r; do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
                 call TYPE3_DET(jh_p,jr_p,jtype3_p)
                call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )   
            EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ prob_s_treat_new_ret(jage, je,jht,jtype3,jtype2, jmar, jMEcat_p, jtype2_p,  jtype3_p, jmar_p ) * (    p_surv(jh_p,je,jmar_p, jage+1) * RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)  +   (1.0d0-p_surv(jh_p,je,jmar_p, jage+1))* (cost_death + Bequest_val(ASP) ))
             end do; end do; end do  ; end do  ; end do 
             call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
	        RET(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,jmar)=UTI+ beta_hat(je)*EXPECTED_RET_PRIME    
       
!Search for optimal 
 
            if (RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar).lt.maxB) exit
			  	    maxB=RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)                  
		        
                
    end do !for jsav for retired ppl
     RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)=maxB 
    
! now find the maximum considering the extremes
    if (RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar).eq.0.0d0) then
        RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)=V_c_all_trpay(1)   
    else 
        RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)=max(RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar), V_c_all_trpay(1))  
    end if 
    
 
end if  ! for getting TR
   	       
RET_PRIME(jhc_ret,jMEcat,jas, jtype1,jh,jr,jtype2,jage,jmar)=RET(jhc_ret,jMEcat,jas,jtype1,jh,jr,jtype2,jmar)

end do ! jas
end do ! hc at age 65
end do !jmar
end do ! jMEcat
end do; end do; end do; end do ;  end do  !loop for   jr, jh, js, jdp, jdu, 	 
end do !jage


end do ! jht
end do ! educ


! we solved only for jfix=1 
RET_PRIME(:,:,:, 3,:,:,:,:,:)=RET_PRIME(:,:,:, 1,:,:,:,:,:)
RET_PRIME(:,:,:, 5,:,:,:,:,:)=RET_PRIME(:,:,:, 1,:,:,:,:,:)

RET_PRIME(:,:,:, 4,:,:,:,:,:)=RET_PRIME(:,:,:, 2,:,:,:,:,:)
RET_PRIME(:,:,:, 6,:,:,:,:,:)=RET_PRIME(:,:,:, 2,:,:,:,:,:)

RET_PRIME(:,:,:, 9,:,:,:,:,:)=RET_PRIME(:,:,:, 7,:,:,:,:,:)
RET_PRIME(:,:,:, 11,:,:,:,:,:)=RET_PRIME(:,:,:, 7,:,:,:,:,:)

RET_PRIME(:,:,:, 10,:,:,:,:,:)=RET_PRIME(:,:,:, 8,:,:,:,:,:)
RET_PRIME(:,:,:, 12,:,:,:,:,:)=RET_PRIME(:,:,:, 8,:,:,:,:,:)

RET_PRIME(:,:,:, 15,:,:,:,:,:)=RET_PRIME(:,:,:, 13,:,:,:,:,:)
RET_PRIME(:,:,:, 17,:,:,:,:,:)=RET_PRIME(:,:,:, 13,:,:,:,:,:)

RET_PRIME(:,:,:, 16,:,:,:,:,:)=RET_PRIME(:,:,:, 14,:,:,:,:,:)
RET_PRIME(:,:,:, 18,:,:,:,:,:)=RET_PRIME(:,:,:, 14,:,:,:,:,:)

!************************************************************************************************************************************
! AGE 64
!************************************************************************************************************************************
! initialize and start all do loops
earnings_interpol=0.0d0 
EXP_L_prime(:,:,:,:,:,:,:,:)=0.0d0
B(:,:,:,:,:,:,:,:)=0.0d0

      
jage=n_ret-1   
PRINT *,'Current Age:',jage     

do je = 1, n_e
do jfix=1, min(n_fix , I_state_fix)
do jht=1, min(n_ht, I_state_ht)
    call TYPE1_DET(je,jht,jfix,jtype1)
    
EXP_L_prime_fix(:,:,:,:,:,:,:)=0.0d0
B_fix(:,:,:,:,:,:,:,:,:,:)=0.0d0
B_1(:,:,:,:,:,:,:,:,:)=0.0d0
B_2(:,:,:,:,:,:,:,:,:)=0.0d0
EXP_L_prime_simple(:,:,:,:,:,:,:)=0.0d0
    
do jhc =1, n_hc_grid     ! from 1 to 5 grid points -- note these correspond to actual years of FT experience from 0 to 39 - given by xfn(jhc,jage)
  
do jtype4=1,n_type4
!find jmar and jos
                if (jtype4.eq.1) then
                jmar=1
                jos=1
            else if (jtype4.eq.2) then
                jmar=2
                jos=1
            else if (jtype4.eq.3) then
                jmar=2
                jos=2
            end if 
do jh=1, n_h ; do jr=1, n_r  
    call TYPE3_DET(jh,jr,jtype3)
do jas=1,n_as ; do jzt=1, n_zt    
do jdu=1,n_du ; do jdp=1,n_dp ; do js=1,n_s
    call TYPE2_DET(jdu,jdp,js,jtype2) !I combine all health shocks into one variable 

do jMEcat=1, n_MEcat  
     
do jo=1, n_o
    
         if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(1)) then 
        jmc=1
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(2)) then
        jmc=2
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(3)) then
        jmc=3
        else
        jmc=4
        end if
        
!initialize
  actual_earn =  0.0d0
  income_bt=0.0d0
  tax_SS=0.0d0
  tax_med=0.0d0
  income_pay(:)=0.0d0
  Total_tax_pay(:)=0.0d0
  base_2_pay(:)=0.0d0
  base_1_pay(:)=0.0d0
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
  V_cf_trpay(:)=0.0d0
  V_c_all_trpay(:)=0.0d0
             SAVE_RET_PRIME_INT(:,:,:,:,:,:)=0.0d0
           SAVE_EXPECTED_B_PRIME_trpay(:,:)=0.0d0
  TR_indic=0
  
!find earnings and income and jhc_ret that determine SS next period

    if (jo.eq.1)  then 
      earnings_interpol=0.0d0
  else 
      call InterpolateScalar ( HK(0 : n_hc), base_earn(0:n_hc,jh,je,jfix,jzt,jo) , xfn(jhc,jage)  , earnings_interpol ) 
  end if 
       ! note base_earn already incorporates the wage penalty for PT , and it is 0 if jo=1
!find earning capacity based on wage offer 
earnings_capacity=0.0d0

 if (xfn(jhc,jage).le.hc_TS(1,je)) then 
    jhc_ret=1
else if (xfn(jhc,jage).le.hc_TS(2,je)) then 
    jhc_ret=2
else 
    jhc_ret=3
end if  
  
actual_earn = max(0.0d0 , ( earnings_interpol * (hrs_offer(jo) -phi(je,jh,jtype2)  )) ) ! earnings net of sick days
income_bt= ((interest-1.0d0))*as(jas,jtype1,jage) +  actual_earn  ! income before tax is equal to interest income plus earnings net of sick days
if ((jo.eq.3).OR.(jo.eq.5)) then ! husband holds insurance -- so deduct from his income only
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) !! own plus spouse's
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else ! subtract from wife. note this can be zero
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) )) !! own plus spouse's
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)))
end if 

     
!determine income group - this matters for marital status probabilities
if (income_bt.le.INC_TS(1)) then 
jinc=1
else if (income_bt.le.INC_TS(2)) then 
jinc=2
else if (income_bt.le.INC_TS(3)) then 
jinc=3
else if (income_bt.le.INC_TS(4)) then 
jinc=4
else 
jinc=5
end if


 !income and taxes assuming you pay. we assume that you always pay this mimimum costs if no shocks
  income_pay(1) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) -tax_SS - tax_med - max(0.d0, oop(jtype2,jh,jage,jMEcat,jo) + oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(1) = max(0.0d0,  (( (income_pay(1)*5200.0d0) - tax_lambda(jmar) *( (income_pay(1)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(1) =  prem_ephi(jo,jtype4) +  oop(jtype2,jh,jage,jMEcat,jo) +  oop_spouse(jage,je,jtype4,jo) +Total_tax_pay(1)  ! total resources subtracted 
  base_1_pay(1) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(1)        ! cash on hand (if you pay the OOP)

!solve for the EV associated with having 0 assets, and getting treated 
        EXPECTED_B_PRIME_C_FLOOR=0.d0  
        RET_PRIME_INT(:,:,:,:,:)=0.0d0
        V_death=0.0d0
        V_death= cost_death + Bequest_val(0.0d0) 

        do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r 
        do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar    
            call TYPE3_DET(jh_p,jr_p,jtype3_p)
        EXPECTED_B_PRIME_C_FLOOR = EXPECTED_B_PRIME_C_FLOOR + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) * Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
        end do; end do                                                                                 
        end do; end do ; end do
        
!solve for the EV associated with having 0 assets, and not getting treated 

     EXPECTED_B_PRIME_C_FLOOR_NT =0.0d0
        do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r 
        do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar   
              call TYPE3_DET(jh_p,jr_p,jtype3_p)
        EXPECTED_B_PRIME_C_FLOOR_NT = EXPECTED_B_PRIME_C_FLOOR_NT + prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p ) * Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)  +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
        end do; end do                                                                                 
        end do; end do ; end do
    
!Calculate value if on cons floor - and get treated  
 if ( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) then  
    CON=cons_floor(jage, jmar,je,jh)
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON,UTI)   ! the utility is that associated with paying and being treated. 
    B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = UTI+ beta_hat(je)*EXPECTED_B_PRIME_C_FLOOR
  
 else 
     !Calculate value if not on cons floor 
!Find max savings
    MAX_SAV_pay(1) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1)  -cons_floor_at(jage, jmar,je,jh) 
    
              jsav_upper_pay(1)=0
          do while ( MAX_SAV_pay(1) .ge.  sav(jo,jsav_upper_pay(1)+1))
           jsav_upper_pay(1)=jsav_upper_pay(1) +1
           if (jsav_upper_pay(1).EQ.n_sav) exit
          end do 
          
               if (jsav_upper_pay(1).lt.1) then
                 PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
               end if 
               
            
      
! Calculate value of consuming the consumption floor and saving all else: ASP= MAX_SAV+as(jas,jtype1,jage) 
     V_cf_trpay(1)=0.0d0
     ASP= MAX_SAV_pay(1) +as(jas,jtype1,jage)
     V_death=0.0d0
     V_death=cost_death + Bequest_val(ASP) 
    call func_uti(jage,1,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)      
    !Find expected value                                           
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
 
       !here, find the lower bound for assets tomorrow that are associated with MAX_SAV_pay(1) +as(jas,jtype1,jage)
      ! the minimum point is "AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1))" but could be higher than this    
      
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1))
          
          
          do while ((MAX_SAV_pay(1) +as(jas,jtype1,jage)) .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
 
      
      
do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
    call TYPE3_DET(jh_p,jr_p,jtype3_p)

    
    
         if (ASP.eq.0.0d0) then
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else if (( m ).eq.n_as) then ! if there is no next point
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                !use the fact we already solved for lower bound
                slope_as=0.0d0
                inter_as=0.0d0
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(m+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((m+1),jtype1,jage+1) -  as((m),jtype1,jage+1)  )
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((m),jtype1,jage+1)
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    
      !***************************
    
    

EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
end do; end do  ; end do;   end do; end do
  V_cf_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(1)

! Calculate value of consuming all and saving nothing
     V_c_all_trpay(1) =0.0d0
     call func_uti(jage,1,je,jh,jo,jtype4,jtype2, (base_1_pay(1)/(1.0d0+taxc)), UTI)  
     V_c_all_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR ! this was solved earlier as the EV of saving nothing and getting treated

     

! BEGIN SAVINGS LOOP  
           maxB_trpay(1)=V_cf_trpay(1) ! initialize to when you save max amount, off the grid
do jsav= jsav_upper_pay(1) , 1, -1 
                 
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
         if (ASP.le.0.d0)  exit ! already did the case where consume all and save 0
    V_death=0.0d0
    V_death= cost_death + Bequest_val(ASP) 
    CON_pay(1)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(1) -sav(jo,jsav) )/(1.0d0+taxc)
    
         if (CON_pay(1).lt.cons_floor(jage, jmar,je,jh)) then 
        PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
        pause
         end if 
    
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON_pay(1),UTI)
!Find expected value                                           
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
    call TYPE3_DET(jh_p,jr_p,jtype3_p)
!  interpolation to get RET_PRIME_INT

         if (( AS_LB(jas,jtype1,jage+1,jo,jsav) ).eq.n_as) then ! if there is no next point
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                !use the fact we already solved for lower bound
                slope_as=0.0d0
                inter_as=0.0d0
            ! slope = (yvect(loc+1)-yvect(loc))/(xvect(loc+1)-xvect(loc))
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jage+1) -  as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)  )
            !inter = yvect(loc) - slope * xvect(loc)
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)
		    !y1 = inter + slope * x1	   
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    
      !***************************
        ! SAVE THIS FOR LATER BECAUSE WE DON'T WANT TO DO THIS TIME CONSUMING LOOP AGAIN WHEN SOLVING FOR OPTIONS TO NOT PAY, NOT TREAT
        SAVE_RET_PRIME_INT(jsav,jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)
        
EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
end do; end do  ; end do;   end do; end do

SAVE_EXPECTED_B_PRIME_trpay(jsav,1)=EXPECTED_B_PRIME_trpay(1) ! we save this EV of getting treated for every jas. we will use it later when solving for value of jtrpay=2 (treat not pay)
     
B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI + beta_hat(je)*EXPECTED_B_PRIME_trpay(1)
           

        if (B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).lt.maxB_trpay(1)) exit
              maxB_trpay(1)=B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)   
                  
                                   
end do   !for jsav  

        B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(1),  V_c_all_trpay(1)) !V_cf_trpay(1),

   
end if   ! for being on cons floor 
 
 !if no shocks, we only need to solve the value if you "pay and get treated"
if (jtype2.eq.1) then      
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
  
 ! *******************************************************************************************************************************************************   
else 
!IF YOU HAVE HEALTH SHOCK, YOU HAVE TO COMPARE not TREAT AND not PAY OPTIONS
!income and taxes if not pay
  
  income_pay(2) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)  -tax_SS - tax_med -  max(0.d0,  oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(2) = max(0.0d0,  (( (income_pay(2)*5200.0d0) - tax_lambda(jmar) *( (income_pay(2)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(2) =  prem_ephi(jo,jtype4) +  oop_spouse(jage,je,jtype4,jo) + Total_tax_pay(2)  ! total resources subtracted 
  base_1_pay(2) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(2)        ! cash on hand (if you do not pay the OOP)
   
!We don't need to consider the option of being on the consumption floor anymore. we already solved for that when we did the first case.
  ! so if that option was ever optimal, it's already been done and captured by B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).
      
!Calculate value if not on cons floor    
 if ( base_1_pay(2) .gt. (cons_floor_at(jage, jmar,je,jh))) then 
 
!Find max savings if you don't pay
   MAX_SAV_pay(2) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(2)  -cons_floor_at(jage, jmar,je,jh) ! if not pay OOP
   
   
             jsav_upper_pay(2)=0
          do while ( MAX_SAV_pay(2) .ge.  sav(jo,jsav_upper_pay(2)+1))
           jsav_upper_pay(2)=jsav_upper_pay(2) +1
           if (jsav_upper_pay(2).EQ.n_sav) exit
          end do 

      
! Calculate value of consuming the consumption floor and saving all else: ASP= MAX_SAV+as(jas,jtype1,jage) 
     V_cf_trpay(2:3)=0.0d0
     ASP= MAX_SAV_pay(2) +as(jas,jtype1,jage)  
     V_death=0.0d0
     V_death=cost_death + Bequest_val(ASP) 
    !Find expected value                                           
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
 
       if (ASP .lt. as(AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2)),jtype1,jage+1)) then 
          pause 
      else  
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2))
          
          do while (ASP .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
      end if 
      
do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
    call TYPE3_DET(jh_p,jr_p,jtype3_p)

  !****************************
         if (ASP.eq.0.0d0) then
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else if (( m ).eq.n_as) then ! if there is no next point
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,m,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                slope_as=0.0d0
                inter_as=0.0d0
            ! slope = (yvect(loc+1)-yvect(loc))/(xvect(loc+1)-xvect(loc))
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(m+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((m+1),jtype1,jage+1) -  as((m),jtype1,jage+1)  )
            !inter = yvect(loc) - slope * xvect(loc)
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(m),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((m),jtype1,jage+1)
		    !y1 = inter + slope * x1	   
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    
      !***************************    
    
! get treated
EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) +  prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)

!not treated
EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) +  prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo)*(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)

end do; end do  ; end do;   end do; end do
   call func_uti(jage,2,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
   V_cf_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(2)
  call func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
   V_cf_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(3)


! Calculate value of consuming all and saving nothing - don't pay, get treated
     V_c_all_trpay(2) =0.0d0
     call func_uti(jage,2,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI)  
     V_c_all_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR ! this was solved earlier as the EV of saving nothing and getting treated
! Calculate value of consuming all and saving nothing - don't pay, do not get treated
     V_c_all_trpay(3) =0.0d0
     V_death=0.0d0
     V_death= cost_death + Bequest_val(0.0d0) 
     call func_uti(jage,3,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI) 
     V_c_all_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 

      
! BEGIN SAVINGS LOOP    
     maxB_trpay(2)=V_cf_trpay(2) ! initiate to when you save max amount
     maxB_trpay(3)=V_cf_trpay(3)
do jsav= jsav_upper_pay(2) , 1, -1 
                    
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
    if (ASP.le.0.d0)  exit ! already did the case where consume all and save 0
    
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
    
    CON_pay(2)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(2) -sav(jo,jsav) )/(1.0d0+taxc)
        
        
         if (CON_pay(2).lt.cons_floor(jage, jmar,je,jh)) then 
          PRINT *,'Error',jas,jtype1,jh,jr,jtype2, jage
          pause
         end if 

                    call  func_uti(jage,2,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(2) = UTI
                    call  func_uti(jage,3,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(3) = UTI     
                                      
!CALCULATE EXPECTED VALUE 
!  we already solved it for some jsav values. only solve it here if we did not solve it before
if ((SAVE_EXPECTED_B_PRIME_trpay(jsav,1)).ne.(0.0d0)) then ! we already solved for the EV of treating. and we can use the saved interpolated values to find value of not treating too.
    B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(2) + beta_hat(je)*SAVE_EXPECTED_B_PRIME_trpay(jsav,1)
  
! solve EV for not pay not treat 
    EXPECTED_B_PRIME_trpay(3)=0.0d0
    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
        call TYPE3_DET(jh_p,jr_p,jtype3_p)   
    EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) + prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* SAVE_RET_PRIME_INT(jsav,jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
    end do; end do  ; end do;   end do; end do
    B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
else ! for jsav we have not considered before, solve EV
 EXPECTED_B_PRIME_trpay(:)=0.d0
 RET_PRIME_INT(:,:,:,:,:)=0.0d0
     do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r;  do jMEcat_p=1, n_MEcat; do jmar_p=1,n_mar
         call TYPE3_DET(jh_p,jr_p,jtype3_p)
     !****************************

         if (( AS_LB(jas,jtype1,jage+1,jo,jsav) ).eq.n_as) then ! if there is no next point
            RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
        else
      
            if (abs(RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) -ASP)<1.d-10) then
                 RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)=RET_PRIME(jhc_ret,jMEcat_p,AS_LB(jas,jtype1,jage+1,jo,jsav),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p)
            else 
                !use the fact we already solved for lower bound
                slope_as=0.0d0
                inter_as=0.0d0
            ! slope = (yvect(loc+1)-yvect(loc))/(xvect(loc+1)-xvect(loc))
		        slope_as =(RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p))   / (  as((AS_LB(jas,jtype1,jage+1,jo,jsav)+1),jtype1,jage+1) -  as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)  )
            !inter = yvect(loc) - slope * xvect(loc)
                inter_as = RET_PRIME(jhc_ret,jMEcat_p,(AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) - slope_as * as((AS_LB(jas,jtype1,jage+1,jo,jsav)),jtype1,jage+1)
		    !y1 = inter + slope * x1	   
                RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) = inter_as + slope_as * ASP	            
            end if
        end if    
      !***************************
 
            EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) + prob_s_treat_new_64(1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death )
            EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) + prob_s_treat_new_64(jmc+1,  je,jht,jtype3,jtype2, jMEcat_p, jtype2_p,  jtype3_p )* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death )
     
     end do; end do  ; end do;   end do; end do
     
           B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(2) + beta_hat(je)*EXPECTED_B_PRIME_trpay(2)
           B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
end if            
 

! if both values are going down, stop and exit
     if ((B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).lt.maxB_trpay(2)).AND.(B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).LT.maxB_trpay(3)))  exit  
         
              if (B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).ge.maxB_trpay(2)) then                              
              maxB_trpay(2)=B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)    
              end if   
              
              if (B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage).ge.maxB_trpay(3)) then                              
              maxB_trpay(3)=B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  ! will save the new max and continue iterating over savings. 
              end if 
              
         
end do   !for jsav    


        B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max( maxB_trpay(2),  V_c_all_trpay(2)) 
        B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max( maxB_trpay(3),  V_c_all_trpay(3)) 

    
        
               
 else ! if on cons floor even when don't pay
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 ! solve if on cons floor and not treat.
    call  func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  =   UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 
 
 ! WHEN WE ALLOW ALL THOSE ON MEDICAID TO GET TREATED, WE HAVE TO CHANGE  B_fix(3,...) TO THE MAX OVER WHAT IT IS NOW, AND THE VALUE OF BEING ON THE CONS FLOOR AND TREATING, BUT ONLY IF YOU QUALIFY FOR IT WHEN YOU PAY.
 if ((I_treat_Mcaid.EQ.1).AND.( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) ) THEN
      B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = max( B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) , B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)) 
 END IF 
 
 
 
 end if   ! for being on cons floor  
end if  ! for jtype2=1 
 !find the max over pay-treat options

!for those with health shocks of type 1 or 2 - so not "can treat and can default"
!among these, the uninsured can't get treated with some probability
! B_1 corresponds to these - I take the max over VF of "treat and pay" and "not treat"
if (jtype2.eq.1) then
B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) )
else ! so here we take prob over can treat or cannot treat (conditional on that you are not "can treat and not pay")
    if ((Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1)).NE.0.0d0) then  
B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) = (Shock_type(2,jo,jh,n_ret-1)/(Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1)))* (max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) ))    + (Shock_type(1,jo,jh,n_ret-1)/(Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1)))* B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1)

    end if 
end if 


!for those who have all 3 options.   
B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1),B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) ,B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) )

end do ! jo
end do ! jMEcat
end do; end do; end do ! for du, dp, s

! Find "EXPECTED_B" over health shocks, ME_cat, and if you have option to not pay and treat.
        EXPECTED_B(:)=0.0d0 
    do jo=1, n_o
            do jtype2=1, n_type2 
            do jMEcat=1, n_MEcat
        
                       
                prob_s=0.0d0
                prob_s= shocks_risk(je,jtype2,jh,jr, jage)*Prob_MEcat(jMEcat) 
   EXPECTED_B(jo)= EXPECTED_B(jo)+  prob_s * ((Shock_type(1,jo,jh,n_ret-1)+Shock_type(2,jo,jh,n_ret-1))* B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1)     + Shock_type(3,jo,jh,n_ret-1) *B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,n_ret-1) )
        
            end do 
            end do 
            B(jo,jtype4,jzt,jas,jtype3,jtype1,jhc,n_ret-1) = EXPECTED_B(jo)
    end do  
  
 ! decide if accept offer and save EXP_L_prime_fix
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,:,jtype4)=EXPECTED_B(1) ! initialize to the non-employed value
    
    If   (EXPECTED_B(5).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,5,jtype4)=EXPECTED_B(5)
    end if
  
    If   (EXPECTED_B(4).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,4,jtype4)=EXPECTED_B(4)
    end if  
    
    If   (EXPECTED_B(3).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,3,jtype4)=EXPECTED_B(3)
    end if
  
    If   (EXPECTED_B(2).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,2,jtype4)=EXPECTED_B(2)
    end if 

end do; end do; end do ; end do   !loop for zt, jas, jr, jh 
end do ! jtype4
end do ! jhc

!Calculate "EXP_L_prime_simple", eliminating jzt and jos, taking probabilities over these.

do jo=1,n_o  ; do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid  
!get back jh because it's needed for probability of wife working
        if ((jtype3.eq.1).OR.(jtype3.eq.2)) then 
        jh=1
        else if ((jtype3.eq.3).OR.(jtype3.eq.4)) then
    jh=2
    else if ((jtype3.eq.5).OR.(jtype3.eq.6)) then
        jh=3
    end if 

    
   ! FOR THOSE NOT EMPLOYED LAST YEAR 
   EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,:,1) =0.0d0
    do jzt=1,n_zt 
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,1) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,1)  + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,1) * p_zt(je,jzt,1)
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,1) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,1)  + p_zt(je,jzt,1) * (EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,2) *prob_os(n_ret-1,je,jh,1,jfix) + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,3) *prob_os(n_ret-1,je,jh,2,jfix))
! for singles, it's the same as jtype=1. for married, take average over the 2 employment states
    end do
   
    ! FOR THOSE EMPLOYED LAST YEAR
       EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,:,2) =0.0d0
    do jzt=1,n_zt 
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,2) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,1,2)  + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,1) * p_zt(je,jzt,2)
EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,2) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,2,2)  + p_zt(je,jzt,2) * (EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,2) *prob_os(n_ret-1,je,jh,1,jfix) + EXP_L_prime_fix(jas,jtype3,jhc,n_ret-1,jzt,jo,3) *prob_os(n_ret-1,je,jh,2,jfix))
    end do
    
end do; end do; end do; end do



do jo=1,n_o  ; do jmar =1,n_mar ; do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid 
EXP_L_prime(jas,jtype1,jtype3,jhc,n_ret-1,jo,jmar,1) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,jmar,1)
EXP_L_prime(jas,jtype1,jtype3,jhc,n_ret-1,jo,jmar,2) =EXP_L_prime_simple(jas,jtype3,jhc,n_ret-1,jo,jmar,2)
end do; end do; end do; end do; end do


end do; end do ; end do  ! jht, jfix, je

 
!******************************************************************************************************************************************
! AGES 25-63  
!******************************************************************************************************************************************

open(14, file='Running Time.txt')

Call cpu_time(time2)
time=(time2-time1)/60
write (14, *) 'Ages 64+ has been running for', time, 'minutes'



do je = 1, n_e
do jfix=1, min(n_fix , I_state_fix)
do jht=1, min(n_ht, I_state_ht)
call TYPE1_DET(je,jht,jfix,jtype1) ! Gives Type1 given educ, fixed productivity and health investment type
B_fix(:,:,:,:,:,:,:,:,:,:)=0.0d0
EXP_L_prime_fix(:,:,:,:,:,:,:)=0.0d0

B_1(:,:,:,:,:,:,:,:,:)=0.0d0
B_2(:,:,:,:,:,:,:,:,:)=0.0d0

ageloop: do jage=n_ret-2, 1, -1   
hkloop: do jhc =1, n_hc_grid  ! do all possible grid points for human capital: 1-5    
    do jh=1, n_h 
do jdu=1,n_du ; do jdp=1,n_dp ; do js=1,n_s 
    call TYPE2_DET(jdu,jdp,js,jtype2)
    
! human capital on  scale 1-100 for tomorrow given my work status jo today 
! this is a real number

jx_p_o(:)=0.0d0
jx_p_o(1) = xfn(jhc,jage)
jx_p_o(2) = min((xfn(jhc,jage)+( max((hrs_pt - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft))) ,  dble(n_hc) )
jx_p_o(3) = jx_p_o(2)
jx_p_o(4) = min( (xfn(jhc,jage)+( max((hrs_ft - phi(je,jh,jtype2)), 0.0d0 ) / (hrs_ft))) ,  dble(n_hc) )
jx_p_o(5) = jx_p_o(4)


do jzhc=1,n_zhc  
  HK_INTERPOL(jzhc,:)  =  max(0.0d0, min( jx_p_o(:) *zhc(jzhc), dble(n_hc) ))  
end do  

!*************************************************************************************************************  
     
do jr=1, n_r 
    call TYPE3_DET(jh,jr,jtype3)  ! Gives Type3 given jh and jr
do jas=1,n_as  
do jtype4=1,n_type4 
!find jmar and jos    
            if (jtype4.eq.1) then
                jmar=1
                jos=1
            else if (jtype4.eq.2) then
                jmar=2
                jos=1
            else if (jtype4.eq.3) then
                jmar=2
                jos=2
            end if    
do jMEcat=1, n_MEcat  
  do jo=1,n_o 
      
      
      ! jo_lagged is 1 if not working, 2 if working. this is used when calculating the expected value tomorrow which depends on the current working status since wages tomorrow depend on lagged employmend. 
      if (jo.eq.1) then 
          jo_lagged=1 
      else 
          jo_lagged=2
      end if 
      
      
     if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(1)) then 
        jmc=1
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(2)) then
        jmc=2
        else if (oop(jtype2,jh,jage,jMEcat,jo).le.MC_TS(3)) then
        jmc=3
        else
        jmc=4
        end if
        

    
do jzt=1, n_zt
    if ((jzt.ne.1).AND.(jo.eq.1)) exit
     
!************************************************************************************************************
!     CALCULATE IF YOU ARE ELIGIBLE FOR TRANSFER AND MAX SAVINGS, AND YOUR INCOME GROUP
!************************************************************************************************************
!initialize
  actual_earn =  0.0d0
  income_bt=0.0d0
  tax_SS=0.0d0
  tax_med=0.0d0
  income_pay(:)=0.0d0
  Total_tax_pay(:)=0.0d0
  base_2_pay(:)=0.0d0
  base_1_pay(:)=0.0d0
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
  V_cf_trpay(:)=0.0d0
  V_c_all_trpay(:)=0.0d0
  m=0.0d0
  n=0.0d0
  SAVE_EXPECTED_B_PRIME_trpay(:,:)=0.0d0 ! this changes with every new state considered. 
SAVE_EXP_L_prime_INT(:,:,:,:,:,:)=0.0d0
  TR_indic=0
  

! earnings and income     
  if (jo.eq.1)  then 
      earnings_interpol=0.0d0
  else 
      call InterpolateScalar ( HK(0 : n_hc), base_earn(0:n_hc,jh,je,jfix,jzt,jo) , xfn(jhc,jage)  , earnings_interpol ) 
  end if 

  
  actual_earn = max(0.0d0 , ( earnings_interpol * (hrs_offer(jo) -phi(je,jh,jtype2) )) ) ! earnings net of sick days
  income_bt= ((interest-1.0d0))*as(jas,jtype1,jage) +  actual_earn  ! income before tax is equal to interest income plus earnings net of sick days
  if ((jo.eq.3).OR.(jo.eq.5)) then ! husband holds insurance -- so deduct from his income only
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else 
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) )) 
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)))
end if 
  
  ! determine income group based on thresholds
if (income_bt.le.INC_TS(1)) then 
jinc=1
else if (income_bt.le.INC_TS(2)) then 
jinc=2
else if (income_bt.le.INC_TS(3)) then 
jinc=3
else if (income_bt.le.INC_TS(4)) then 
jinc=4
else 
jinc=5
end if

!income and taxes assuming you pay. 
  income_pay(1) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4) -tax_SS - tax_med - max(0.d0, oop(jtype2,jh,jage,jMEcat,jo) + oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(1) = max(0.0d0,  (( (income_pay(1)*5200.0d0) - tax_lambda(jmar) *( (income_pay(1)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 ))  +  tax_SS + tax_med  ! tax 
  !Tax_inc(jmar,income_pay(1)) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(1) =  prem_ephi(jo,jtype4) +  oop(jtype2,jh,jage,jMEcat,jo) +  oop_spouse(jage,je,jtype4,jo) +Total_tax_pay(1)  ! total resources subtracted 
  base_1_pay(1) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(1)        ! cash on hand (if you pay the OOP)
  
! solve for the EV associated with having 0 assets, and getting treated, and also 0 assets but not getting treated.
  EXPECTED_B_PRIME_C_FLOOR=0.0d0
  EXPECTED_B_PRIME_C_FLOOR_NT =0.0d0 ! if you don't get treated
  EXP_L_prime_INT=0.0d0
  check=0.0d0

        do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
               call TYPE3_DET(jh_p,jr_p,jtype3_p) 
            !INTERPOLATE FOR HUMAN CAPITAL
            if (HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1) .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(1,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(1,jtype3_p, (HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(1,jtype3_p, HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1), jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1))+1, jage+1)  -   xfn((HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)), jage+1))
		    inter_hc = EXP_L_prime_simple(1,jtype3_p, HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1), jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            

            EXPECTED_B_PRIME_C_FLOOR= EXPECTED_B_PRIME_C_FLOOR      + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) * (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * (cost_death + Bequest_val(0.0d0) )  ) 
            EXPECTED_B_PRIME_C_FLOOR_NT=EXPECTED_B_PRIME_C_FLOOR_NT + prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) * (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * (cost_death + Bequest_val(0.0d0) )  ) 
        end do ; end do ; end do; end do ; end do    

        
                 
!Calculate value if on cons floor - here you get treated  
if ( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) then
 
    CON=cons_floor(jage, jmar,je,jh)
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON,UTI)   ! the utility is that associated with paying and being treated. 
    B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = UTI+ beta_hat(je)*EXPECTED_B_PRIME_C_FLOOR
           

else   
     !Calculate value if not on cons floor     
!Find max savings
    MAX_SAV_pay(1) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1)  -cons_floor_at(jage, jmar,je,jh) 

                !find lower bound 
          jsav_upper_pay(1)=0
          do while ( MAX_SAV_pay(1) .ge.  sav(jo,jsav_upper_pay(1)+1))
           jsav_upper_pay(1)=jsav_upper_pay(1) +1
           if (jsav_upper_pay(1).EQ.n_sav) exit
          end do 
 
    
 ! Calculate value of consuming the consumption floor and saving all else: ASP= MAX_SAV+as(jas,jtype1,jage) 
     V_cf_trpay(1)=0.0d0
     ASP= MAX_SAV_pay(1) +as(jas,jtype1,jage)
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
    call func_uti(jage,1,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)      
    !Find expected value                                           
 EXPECTED_B_PRIME_trpay(:)=0.d0   
 EXP_L_prime_INT=0.0d0
     
       !here, find the lower bound for assets tomorrow that are associated with MAX_SAV_pay(1) +as(jas,jtype1,jage)
      ! the minimum point is "AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1))" but could be higher than this    
     
      if (ASP .lt. as(AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1)),jtype1,jage+1)) then 
          pause 
      else  
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1))
          
          do while (ASP .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
      end if 
      
     do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)  
          
            
                !***********************************
                           
! Interpolate for assets and human capital to get "EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)"
        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)

        if (as(n_as,jtype1,jage+1) <= ASP) then
            !INTERPOLATE FOR HUMAN CAPITAL
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn(n, jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 78 
 
        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
            !INTERPOLATE FOR ASSETS. 
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 78             
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 78                
  
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) + &
                mx1*(1-mx2)*yt(2) + &
                (1-mx1)*mx2*yt(3) + &
                mx1*mx2*yt(4) 

78    EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo)  *(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death )
    end do; end do  ; end do;   end do; end do
      V_cf_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(1)

! Calculate value of consuming all and saving nothing
     V_c_all_trpay(1) =0.0d0
     call func_uti(jage,1,je,jh,jo,jtype4,jtype2, base_1_pay(1)/(1.0d0+taxc), UTI)  
     V_c_all_trpay(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR ! this was solved earlier as the EV of saving nothing and getting treated    
    


! BEGIN SAVINGS LOOP 
 maxB_trpay(1)=V_cf_trpay(1)
do jsav= jsav_upper_pay(1) , 1, -1 
                 
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
        if (ASP.le.0.d0)  exit ! we already did the case where consume all and save 0
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
    CON_pay(1)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(1) -sav(jo,jsav) )/(1.0d0+taxc)
    
    
    call  func_uti(jage,1,je,jh,jo,jtype4,jtype2,CON_pay(1),UTI)

!Find expected value
 EXPECTED_B_PRIME_trpay(:)=0.d0   
 EXP_L_prime_INT=0.0d0  
 VF=0.0d0
      do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)                    
          
! Interpolate for assets and human capital to get "EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)"
        m=AS_LB(jas,jtype1,jage+1,jo,jsav)
        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1) 
        

        if (as(n_as,jtype1,jage+1) <= ASP) then
            !INTERPOLATE FOR HUMAN CAPITAL
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn(n, jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 79 

        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
            !INTERPOLATE FOR ASSETS. 
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 79             
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
             !INTERPOLATE FOR ASSETS.
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 79                
  
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) + &
                mx1*(1-mx2)*yt(2) + &
                (1-mx1)*mx2*yt(3) + &
                mx1*mx2*yt(4) 
     
     !  *******************************************
       
        ! SAVE THIS FOR LATER BECAUSE WE DON'T WANT TO DO THIS TIME CONSUMING LOOP AGAIN WHEN SOLVING FOR OPTIONS TO NOT PAY, NOT TREAT
        ! WE ONLY  NEED TO ADD jas TO THE LIST OF VARIABLES IN EXP_L_prime_INT. This is because ASP only depends on jas. 
        
79        SAVE_EXP_L_prime_INT(jsav,jh_p,jr_p, jo_p,jmar_p,jzhc_p) =  EXP_L_prime_INT
  	    	        

            EXPECTED_B_PRIME_trpay(1)= EXPECTED_B_PRIME_trpay(1) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)    
      end do ; end do ; end do   ; end do ; end do 

SAVE_EXPECTED_B_PRIME_trpay(jsav,1)=EXPECTED_B_PRIME_trpay(1) ! we save this EV of getting treated for every jas. we will use it later when solving for value of jtrpay=2 (treat not pay)
     
VF=UTI + beta_hat(je)*EXPECTED_B_PRIME_trpay(1)


     if (VF.lt.maxB_trpay(1)) exit
        
        if ((VF).ge.(maxB_trpay(1))) then                              
              maxB_trpay(1)=VF  
        end if  

                                
end do  !jsav 
    
! now find the maximum considering the extremes 
 
        B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(1),  V_c_all_trpay(1)) 

   
end if   ! for being on cons floor     
    
 !if no shocks, we only need to solve the value if you "pay and get treated"
if (jtype2.eq.1) then      
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
  
 ! *******************************************************************************************************************************************************   
else 
!IF YOU HAVE HEALTH SHOCK, YOU HAVE TO COMPARE not TREAT AND not PAY OPTIONS
!income and taxes if not pay
    
  income_pay(2) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo,jtype4)  -tax_SS - tax_med -  max(0.d0,  oop_spouse(jage,je,jtype4,jo) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(2) = max(0.0d0,  (( (income_pay(2)*5200.0d0) - tax_lambda(jmar) *( (income_pay(2)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 ))  +  tax_SS + tax_med  ! tax 
  base_2_pay(2) =  prem_ephi(jo,jtype4) +  oop_spouse(jage,je,jtype4,jo) + Total_tax_pay(2)  ! total resources subtracted 
  base_1_pay(2) =  interest*as(jas,jtype1,jage) + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(2)        ! cash on hand (if you do not pay the OOP)

!We don't need to consider the option of being on the consumption floor anymore. we already solved for that when we did the first case, jtrpay. 
      
!Calculate value if not on cons floor    
 if ( base_1_pay(2) .gt. (cons_floor_at(jage, jmar,je,jh))) then 
 
!Find max savings if you don't pay
   MAX_SAV_pay(2) = income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(2)  -cons_floor_at(jage, jmar,je,jh) ! if not pay OOP
        
          !find lower bound 
          jsav_upper_pay(2)=0
          do while ( MAX_SAV_pay(2) .ge.  sav(jo,jsav_upper_pay(2)+1))
           jsav_upper_pay(2)=jsav_upper_pay(2) +1
           if (jsav_upper_pay(2).EQ.n_sav) exit
          end do 
          
      
! Calculate value of consuming the consumption floor and saving all else: ASP= MAX_SAV+as(jas,jtype1,jage) - here, you don't pay oop and get treated; and also solve when  you don't get treated.
     V_cf_trpay(2)=0.0d0
     V_cf_trpay(3)=0.0d0
     ASP= MAX_SAV_pay(2) +as(jas,jtype1,jage)  
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
     !Find expected value                                           
 EXPECTED_B_PRIME_trpay(:)=0.d0   
 EXP_L_prime_INT=0.0d0
 
        !here, find the lower bound for assets tomorrow that are associated with MAX_SAV_pay(1) +as(jas,jtype1,jage)
      ! the minimum point is "AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(1))" but could be higher than this    
     
      if (ASP .lt. as(AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2)),jtype1,jage+1)) then 
          pause 
      else  
          m=AS_LB(jas,jtype1,jage+1,jo,jsav_upper_pay(2))
          
          do while (ASP .ge. as((m+1),jtype1,jage+1))
           m=m +1
           if (m.EQ.n_as) exit
           end do 
      end if 
      
 
     do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)  
    !***********************************
                  
          
! Interpolate for assets and human capital to get "EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)"
        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)
 
        
        if (ASP ==0.0d0 .and. HK_INTERPOL(jzhc_p,jo) == 0.0d0) then
           EXP_L_prime_INT= EXP_L_prime_simple(1,jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)
           go to 87
        end if
        if (ASP ==0.0d0) then	
                        !INTERPOLATE FOR HUMAN CAPITAL
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(1,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(1,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(1,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(1,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if      
            go to 87             
        end if
        if (as(n_as,jtype1,jage+1) <= ASP) then
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 87 

        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
            !INTERPOLATE FOR ASSETS. 
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 87             
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
             !INTERPOLATE FOR ASSETS. 
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 87                
 
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) + &
                mx1*(1-mx2)*yt(2) + &
                (1-mx1)*mx2*yt(3) + &
                mx1*mx2*yt(4) 

     !  *******************************************          

            
87    EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
    EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) + prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo)*(p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)

     end do; end do  ; end do;   end do; end do
         call func_uti(jage,2,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_cf_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(2)  
         call func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_cf_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_trpay(3)
  
! Calculate value of consuming all and saving nothing - don't pay, get treated; and also for don't get treated
     V_c_all_trpay(2) =0.0d0
     call func_uti(jage,2,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI)  
     V_c_all_trpay(2) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR ! this was solved earlier as the EV of saving nothing and getting treated  
  
  !if not get treated
       V_c_all_trpay(3) =0.0d0
     call func_uti(jage,3,je,jh,jo,jtype4,jtype2, base_1_pay(2)/(1.0d0+taxc), UTI) 
  V_c_all_trpay(3) = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 
  
   
! BEGIN SAVINGS LOOP   
  !we initiate the max to be the value associated with saving all you can save, ie, V_cf_trpay()
  maxB_trpay(2)=V_cf_trpay(2)
  maxB_trpay(3)=V_cf_trpay(3)
do jsav= jsav_upper_pay(2) , 1, -1 
                    
    ASP=  as(jas,jtype1,jage)+ sav(jo,jsav)
     V_death=0.0d0
     V_death= cost_death + Bequest_val(ASP) 
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
    CON_pay(2)= ( income_bt  + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(2) -sav(jo,jsav) )/(1.0d0+taxc)
     

                    call  func_uti(jage,2,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(2) = UTI
                    call  func_uti(jage,3,je,jh,jo,jtype4,jtype2,CON_pay(2),UTI)
                UTI_trpay(3) = UTI   
                
 VF_tr(2:3)=0.0d0                                        
!CALCULATE EXPECTED VALUE 
! we already solved it for some jsav values.  only solve it here if we did not solve it before
if ((SAVE_EXPECTED_B_PRIME_trpay(jsav,1)).ne.(0.0d0)) then ! we already solved for the EV of treating. and we can use the saved interpolated values to find value of not treating too.
    VF_tr(2)=UTI_trpay(2) + beta_hat(je)*SAVE_EXPECTED_B_PRIME_trpay(jsav,1)
  
! solve EV for not pay not treat 
    EXPECTED_B_PRIME_trpay(3)=0.0d0
       do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)    
        EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) +  prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *(p_surv(jh_p,je,jmar_p, jage+1)* SAVE_EXP_L_prime_INT(jsav,jh_p,jr_p, jo_p,jmar_p,jzhc_p) +  (1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death)
        end do; end do  ; end do;   end do; end do    
            
    VF_tr(3)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
                
                
                
else 
 ! for jsav we have not considered before, solve EV               
            EXP_L_prime_INT=0.0d0
            EXPECTED_B_PRIME_trpay(:)=0.0d0   
             m=AS_LB(jas,jtype1,jage+1,jo,jsav)
             
     do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar  ;  do jzhc_p=1,n_zhc
            call TYPE3_DET(jh_p,jr_p,jtype3_p)                    
          
! Interpolate for assets and human capital to get "EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)"
        
        n=HK_LB(jhc,je,jtype2,jh,jzhc_p,jo,jage+1)
        
        if (as(n_as,jtype1,jage+1) <= ASP) then
            !INTERPOLATE FOR HUMAN CAPITAL
            if (n .ge. n_hc_grid) then
                EXP_L_prime_INT=EXP_L_prime_simple(n_as,jtype3_p, n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)
            else 
                slope_hc =0.0d0
                inter_hc =0.0d0
            slope_hc = (EXP_L_prime_simple(n_as,jtype3_p, (n+1), jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged))  /( xfn((n)+1, jage+1)  -   xfn((n), jage+1))
		    inter_hc = EXP_L_prime_simple(n_as,jtype3_p, n, jage+1,jo_p,jmar_p,jo_lagged)  - slope_hc *  xfn((n), jage+1)
		    EXP_L_prime_INT = inter_hc + slope_hc * HK_INTERPOL(jzhc_p,jo)
            end if 
            go to 88 
        end if
        
         
        if (HK_INTERPOL(jzhc_p,jo) == 0.0d0) then	
            !INTERPOLATE FOR ASSETS. 
                slope_as=0.0d0
                inter_as=0.0d0
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged))/( as((m+1),jtype1,jage+1)-as(m,jtype1,jage+1)  )
		    inter_as = EXP_L_prime_simple((m),jtype3_p,1, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as(m,jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP		
            go to 88            
        end if        
        if (xfn(n_hc_grid, jage+1) <= HK_INTERPOL(jzhc_p,jo)) then
                slope_as=0.0d0
                inter_as=0.0d0
             !INTERPOLATE FOR ASSETS. 
            slope_as = (EXP_L_prime_simple((m+1),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged)-EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged))/(as(m+1,jtype1,jage+1)-as((m),jtype1,jage+1))
		    inter_as = EXP_L_prime_simple((m),jtype3_p,n_hc_grid, jage+1,jo_p,jmar_p,jo_lagged) - slope_as * as((m),jtype1,jage+1)
		    EXP_L_prime_INT = inter_as + slope_as * ASP	
            go to 88                
        end if        
      

	    yt(1) = EXP_L_prime_simple(m,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(2) = EXP_L_prime_simple(m+1,jtype3_p,n, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(3) = EXP_L_prime_simple(m,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
	    yt(4) = EXP_L_prime_simple(m+1,jtype3_p,n+1, jage+1,jo_p,jmar_p,jo_lagged)
        
        if (abs(ASP-as(m,jtype1,jage+1)) < 1.d-10) then
            mx1 = 0.0d0
        else
	        mx1 = (ASP-as(m,jtype1,jage+1))/(as(m+1,jtype1,jage+1)-as(m,jtype1,jage+1))	
	    end if
	    if (abs(HK_INTERPOL(jzhc_p,jo) - xfn(n,jage+1)) < 1.d-10) then
	        mx2 = 0.0d0
	    else
	        mx2 = (HK_INTERPOL(jzhc_p,jo) -xfn(n,jage+1))/(xfn(n+1,jage+1)-xfn(n,jage+1))
        end if
        	
	    EXP_L_prime_INT =    (1-mx1)*(1-mx2)*yt(1) +  mx1*(1-mx2)*yt(2) + (1-mx1)*mx2*yt(3) +  mx1*mx2*yt(4) 

     !  *******************************************
 	    	            
            
88            prob_d=0.0d0
            prob_d=(1.0d0- p_surv(jh_p,je,jmar_p, jage+1)) * V_death
            
            EXPECTED_B_PRIME_trpay(2)= EXPECTED_B_PRIME_trpay(2) + prob_s_treat_new(1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p) *Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  prob_d)
            EXPECTED_B_PRIME_trpay(3)= EXPECTED_B_PRIME_trpay(3) +  prob_s_treat_new(jmc+1,jage, jo,je,jht,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo) *  (p_surv(jh_p,je,jmar_p, jage+1)* EXP_L_prime_INT +  prob_d)

     end do ; end do ; end do   ; end do ; end do 
     

           VF_tr(2)=UTI_trpay(2) + beta_hat(je)*EXPECTED_B_PRIME_trpay(2)
           VF_tr(3)=UTI_trpay(3) + beta_hat(je)*EXPECTED_B_PRIME_trpay(3)
           
end if 



! if both values are going down, stop and exit
     if (((VF_tr(2)).lt.(maxB_trpay(2))).AND.(VF_tr(3).LT.maxB_trpay(3)))  exit  
         
              if (VF_tr(2).ge.maxB_trpay(2)) then                              
              maxB_trpay(2)=VF_tr(2)   
              end if   
              
              if (VF_tr(3).ge.maxB_trpay(3)) then                              
              maxB_trpay(3)=VF_tr(3)  ! will save the new max and continue iterating over savings. 
              end if 
                      
            
         
end do   !for jsav    

   
    
    ! now find the maximum considering the extremes 
  
        B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(2),  V_c_all_trpay(2)) 
        B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)=max(maxB_trpay(3),  V_c_all_trpay(3)) 

    
               
 
 else ! if I am on the consumption floor even when i don't pay, then assign value equal to just that of paying and being on the cons floor. we need this bc we'll take the max over the 3 values, so we need them all to not be 0.
 B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) 
 !  solve value of being on cons floor and not treating
          call func_uti(jage,3,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
 B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  = UTI + beta_hat(je)* EXPECTED_B_PRIME_C_FLOOR_NT 

          
 end if   ! for not being on cons floor
 
 if ((I_treat_Mcaid.EQ.1).AND.( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) ) THEN
      B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = max( B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) , B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)) 
 END IF 
 
end if  ! for jtype2=1  


if (jtype2.eq.1) then
 B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) ) 
else 
    if ((Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage)).NE.0.0d0) then 
B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) = (Shock_type(2,jo,jh,jage)/(Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage)))* max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage), B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) ) + (Shock_type(1,jo,jh,jage)/(Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage))) *  B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)
end if
end if 

  
 

B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) =  max( B_fix(1,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage), B_fix(2,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) ,B_fix(3,jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) )
  
end do ! jzt
end do ! jo
B_1(1,jMEcat,jtype4,2:n_zt,jas,jtype3,jtype2,jhc,jage) = B_1(1,jMEcat,jtype4,1,jas,jtype3,jtype2,jhc,jage) !for those who don't get an offer, zt does not matter, so set all equal to when zt =1
B_2(1,jMEcat,jtype4,2:n_zt,jas,jtype3,jtype2,jhc,jage) = B_2(1,jMEcat,jtype4,1,jas,jtype3,jtype2,jhc,jage) 

end do ! jMEcat
end do ! jmar
end do ! jas
end do ! jr
end do; end do; end do ! for  du, dp, s


do jas =1, n_as ; do jr=1, n_r ;  do jzt=1, n_zt ;  do jtype4=1, n_type4
 call TYPE3_DET(jh,jr,jtype3)  ! Gives Type3 given jh and jr

! Find "EXPECTED_B" over health shocks, ME_cat
        EXPECTED_B(:)=0.0d0 
    do jo=1, n_o
            do jtype2=1, n_type2 
            do jMEcat=1, n_MEcat
                       
                prob_s=0.0d0
                prob_s= shocks_risk(je,jtype2,jh,jr, jage)*Prob_MEcat(jMEcat)
   EXPECTED_B(jo)= EXPECTED_B(jo)+  prob_s * (   (Shock_type(1,jo,jh,jage)+Shock_type(2,jo,jh,jage)) * B_1(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage) + Shock_type(3,jo,jh,jage) * B_2(jo,jMEcat,jtype4,jzt,jas,jtype3,jtype2,jhc,jage)  )
            end do 
            end do 
            
            B(jo,jtype4,jzt,jas,jtype3,jtype1,jhc,jage) = EXPECTED_B(jo) ! this is all we need to save. the expected value given the state at the beginning of the period.

    
    end do  ! jo
              
           
 ! decide if accept offer and save EXP_L_prime_fix
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,:,jtype4)=EXPECTED_B(1) ! initialize to the non-employed value
    
    If   (EXPECTED_B(5).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,5,jtype4)=EXPECTED_B(5)
    end if
  
    If   (EXPECTED_B(4).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,4,jtype4)=EXPECTED_B(4)
    end if  
    
    If   (EXPECTED_B(3).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,3,jtype4)=EXPECTED_B(3)
    end if
  
    If   (EXPECTED_B(2).gt.EXPECTED_B(1)) then  
    EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,2,jtype4)=EXPECTED_B(2)
    end if 

    
end do ; end do; end do ; end do   !loop for  jas, jr, jzt  , jtype4	
end do ! jh
end do hkloop


!Get "EXP_L_prime_simple" which does not depend on jos (wife working) and jzt. We take probabilities over these. 
do jo=1,n_o  ;  do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid   
! get back jh here because it's needed for probability of wife working
        if ((jtype3.eq.1).OR.(jtype3.eq.2)) then 
        jh=1
        else if ((jtype3.eq.3).OR.(jtype3.eq.4)) then
    jh=2
    else if ((jtype3.eq.5).OR.(jtype3.eq.6)) then
        jh=3
    end if 
    
    ! NOT EMPLOYED LAST YEAR
   EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,:,1) =0.0d0
    do jzt=1,n_zt 
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,1) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,1)  + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,1) * p_zt(je,jzt,1)
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,1) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,1)  +   p_zt(je,jzt,1)* (EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,2)*prob_os(jage,je,jh,1,jfix)   + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,3)*prob_os(jage,je,jh,2,jfix) ) 
    end do
   
    !EMPLOYED LAST YEAR
       EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,:,2) =0.0d0
    do jzt=1,n_zt 
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,2) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,1,2)  + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,1) * p_zt(je,jzt,2)
EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,2) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,2,2)  +   p_zt(je,jzt,2)* (EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,2)*prob_os(jage,je,jh,1,jfix)   + EXP_L_prime_fix(jas,jtype3,jhc,jage,jzt,jo,3)*prob_os(jage,je,jh,2,jfix) ) 
    end do
    
end do; end do; end do; end do 

end do ageloop 



do jo=1,n_o  ; do jmar =1,n_mar ; do jas=1,n_as ; do jtype3 =1,n_type3 ;  do jhc =1,n_hc_grid  ; do jage=1,n_ret-2
EXP_L_prime(jas,jtype1,jtype3,jhc,jage,jo,jmar,1) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,jmar,1)
EXP_L_prime(jas,jtype1,jtype3,jhc,jage,jo,jmar,2) =EXP_L_prime_simple(jas,jtype3,jhc,jage,jo,jmar,2)
end do; end do; end do; end do; end do; end do


end do ! jht
end do ! jfix
end do  ! je

Call cpu_time(time4)
time=(time4-time1)/60
write (14, *) 'Solving VF has been running for', time, 'minutes'



!***********************************************************************
! SIMULATE ECONOMY,
 Call SIMULATION(EXP_L_prime, B,  RET_PRIME ,  jExperimentNumber )
!**************************************************************

end do ! for experiments


Call cpu_time(time5)
time=(time5-time4)/60
write (14, *) 'Simulation has been running for', time, 'minutes'
time =(time5-time1)/60
write (14, *) 'This program has been running for', time, 'minutes'
close(14)
End Program VALUEFUNCTION